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1  Executive  Summary 


Practical  scattering  problems  frequently  involve  thin  metal  objects  with 
open  edges  that  can  be  modeled  as  perfectly  conducting  open  surface  tar¬ 
gets.  Examples  include  airborne  vehicles  with  thin  fins,  horn  antennas,  and 
high  frequency  circuits  comprising  thin  strips  of  metallization  on  multi-layer 
substrates.  The  standard  integral  equation  formulations  for  obtaining  fre¬ 
quency  domain  solutions  to  such  open  surface  scattering  problems  often  lead 
to  linear  systems  of  equations  that  are  impractical  to  solve  numerically  -  the 
linear  system  is  either  too  large  to  solve  in  a  reasonable  amount  of  time  by 
direct  inversion  or  too  ill  conditioned  to  solve  by  iterative  means. 

The  most  commonly  used  integral  equation  for  solving  open  surface  scat¬ 
tering  problems  is  the  electric  field  integral  equation  (EFIE).  It  is  the  only 
integral  equation  that  can  be  applied  to  both  open  and  closed  targets.  The 
only  practical  way  to  solve  large  scattering  problems  is  using  iterative  solvers 
in  conjunction  with  operator  decomposition  methods  such  as  the  fast  mul¬ 
tipole  method  (FMM).  Unfortunately,  the  EFIE  is  horribly  ill  conditioned 
so  the  iteration  count  is  uncontrollable  and  can  be  very  large,  thus  negating 
any  advantage  that  may  accrue  from  use  of  a  fast  method. 

The  objective  of  this  program  was  to  produce  a  new  integral  equation 
that  is  inherently  well  conditioned  and  suitable  for  solving  open  surface  scat¬ 
tering  problems.  We  achieved  this  goal  by  developing  and  implementing  an 
analytic  preconditioner  for  the  EFIE.  We  tested  the  efficacy  of  the  new  for¬ 
mulation  on  simple  test  problems,  finding  that  it  stabilized  and  dramatically 
improved  iterative  solver  performance  as  compared  to  conventional  numer¬ 
ical  preconditioning  methods.  For  example,  on  the  EletroMagnetic  Code 
Consortium  (EMCC)  triangle-circle  target  the  analytically  preconditioned 
EFIE  only  required  ~30  iterations  to  achieve  convergence  as  compared  to 
~50000  for  the  standard  EFIE  with  a  conventional  block  diagonal  precon¬ 
ditioner.  In  stark  contrast  to  the  EFIE,  the  eigenvalue  spectrum  of  the  new 
formulation  has  a  distinctly  second  kind  character.  Consequently,  the  com¬ 
puted  source  distribution  (i.e. ,  equivalent  surface  current)  is  also  noticeably 
smoother. 

Our  original  plan  for  realizing  a  numerical  implementation  of  the  analyt¬ 
ically  preconditioned  EFIE  involved  using  a  2d  extension  to  the  Id  Poincare- 
Bertrand  identity  (PBI)  in  order  to  make  the  equation  explicitly  second  kind. 
In  the  process  of  numerically  validating  the  2d  PBI,  however,  we  found  a 
simpler  algorithm  that  accurately  evaluates  the  required  double  integral  op¬ 
erators  in  much  less  time.  In  contrast  to  previously  implemented  methods, 
the  spectrum  of  the  discretized  operator  produced  by  this  method  faithfully 
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reproduces  that  of  the  continuous  operator  even  for  high  spatial  frequency 
eigennrodes. 

One  of  the  surprising  outcomes  of  this  program  is  the  observation  that, 
for  certain  open  surface  PEC  targets,  the  spectrum  of  the  analytically  pre¬ 
conditioned  EFIE  (APEFIE)  includes  one  or  two  zero  eigenvalues.  We  en¬ 
deavored  to  analytically  derive  the  spectrum  of  the  APEFIE  on  a  PEC  disk 
to  determine  the  source  of  the  zero  eigenvalues,  but  were  unable  to  obtain 
definitive  results  before  the  conclusion  of  the  program.  Nonetheless,  our 
report  contains  a  summary  of  the  analytical  results  we  did  obtain  for  open 
surface  scattering  from  canonical  targets  in  the  2d  and  3d  scalar /acoustic 
cases  and  the  3d  vector/EM  case. 

2  Introduction 

The  electric  field  integral  equation  (EFIE)  is  an  essential  element  in  the 
arsenal  of  tools  for  doing  frequency  domain  electromagnetic  scattering  cal¬ 
culations.  It  is  the  only  integral  equation  that  can  be  applied  to  both  open 
and  closed  surface  PEC  targets 

The  EFIE  has  a  serious  problem,  however  -  it  is  notoriously  ill  con¬ 
ditioned.  When  the  system  of  linear  equations  generated  by  the  EFIE  is 
solved  by  an  iterative  method,  the  iteration  count  cannot  be  controlled. 
This  compromises  the  effectiveness  of  acceleration  methods  such  as  the  fast 
multipole  method  (FMM)  because  such  methods  necessarily  rely  on  iterative 
solvers.  Even  when  a  direct  solver  (e.g.,  LU  decomposition)  can  be  used,  the 
combination  of  a  seriously  ill  conditioned  linear  system  and  finite  precision 
computer  arithmetic  can  degrade  the  accuracy  of  the  computed  solution. 

The  goal  of  this  program  is  to  find  an  analytic  preconditioner  for  the 
EFIE  that 

•  converts  the  EFIE  to  a  second  kind  integral  operator, 

•  applies  to  both  open  and  closed  surface  PEC  targets,  and 

•  is  amenable  to  efficient  numerical  implementation. 

Our  approach  is  based  on  two  observations. 

The  first  is  the  EFIE  preconditioning  method  described  in  [1],  When 
applied  to  closed  surface  scattering  problems,  this  analytic  preconditioner 
is  guaranteed  to  convert  the  EFIE  to  a  second  kind  integral  equation.  The 
computational  cost  of  this  approach  is  predictable  and  reasonable.  However, 
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our  original  implementation  of  this  approach  did  not  produce  a  second  kind 
integral  equation  when  applied  to  an  open  surface  scatterer. 

The  second  observation  is  a  2d  extension  [2]  [3]  of  a  well  established,  but 
not  so  widely  known,  Id  integral  identity  called  the  Poincare-Bertrand  iden¬ 
tity  (PBI)  [4].  Analytic  preconditioners  for  the  2d  and  3d  scalar  analogs  of 
the  EFIE  produce  double  integral  operators  that  exhibit  second  kind  behav¬ 
ior.  However,  the  second  kind  behavior  may  be  lost  as  a  result  of  numerical 
discretization.  The  Id  and  2d  Poincare-Bertrand  identities  provide  a  way 
to  write  these  double  integral  operators  in  an  explicitly  second  kind  form, 

i.e. ,  as  the  sum  of  a  constant  term  and  a  compact  operator,  which  makes 
it  trivial  to  retain  second  kind  behavior  in  the  numerical  implementation. 
Furthermore,  the  Poincare-Bertrand  identities  can  be  used  both  with  open 
and  closed  surface  formulations. 

Work  on  the  program  was  divided  into  three  tasks: 

1.  Numerically  validate  the  2d  PBI  on  simple  2d  surfaces. 

2.  Use  the  2d  PBI  to  realize  an  explicitly  second  kind  analytically  pre¬ 
conditioned  EFIE  (APEFIE)  for  closed  surface  PEC  targets. 

3.  Use  the  2d  PBI  to  realize  an  explicitly  second  kind  APEFIE  for  open 
surface  PEC  targets. 

We  completed  all  three  tasks  and  achieved  the  stated  objectives.  In  the 
course  of  our  work,  however,  we  found  a  faster,  more  robust  alternative  to 
the  2d  PBI  for  achieving  second  kind  behavior  in  a  numerical  realization  of 
the  APEFIE.  The  alternative  method  involves  overdiscretizing  the  double 
integral  evaluations  at  intermediate  sample  points.  An  intermediate  sample 
point  density  about  double  that  of  the  field/source  point  discretization  was 
generally  sufficient  to  avoid  compromising  accuracy.  We  used  this  alternative 
method  in  our  numerical  demonstrations. 

Notes: 

•  EFIE  ill  conditioning  has  a  local  cause  and  a  global  cause.  This  pro¬ 
gram  is  only  concerned  with  addressing  the  local  cause.  It  does  not 
directly  address  resonances,  whether  spurious  or  physical,  which  fall 
into  the  global  cause  category. 

•  The  primary  motivation  for  developing  a  well  conditioned  counterpart 
to  the  EFIE  is  to  facilitate  the  use  of  fast  methods  such  as  the  FMM 
for  solving  large  scattering  problems  involving  open  PEC  surfaces. 


3 


The  goal  of  this  program  was  limited  to  finding  a  suitable  integral 
equation  formulation  and  demonstrating  its  potential  on  small  scale, 
open  surface  scattering  problems.  Since  the  underlying  cause  of  the 
problem  is  a  local  one,  the  methods  reported  here  should  be  just  as 
effective  on  large  scale  scattering  problems  as  they  were  for  the  small 
test  targets  investigated  here. 

We  believe  that  this  new  formulation  can  have  far-reaching  impact  be¬ 
cause  it  removes  a  major  impediment  to  the  practical  use  of  fast  methods 
for  solving  a  wide  variety  of  important  scattering  problems.  The  types 
of  problems  that  would  be  affected  range  from  predicting  the  radar  cross 
section  (RCS)  of  airborne  vehicles  with  thin  fins  or  antennas  to  modeling 
microwave  circuits  comprising  thin  strips  of  metallization  on  multi-layer  sub¬ 
strates.  Note  that  the  new  formulation  also  applies  to  3d  scalar  scattering  so 
it  can  be  used  to  improve  solution  methods  for  acoustic  scattering  problems 
as  well. 


3  Poincare-Bertrand  Identity 

The  Poincare-Bertrand  identities  (PBI)  [3]  [4]  provide  a  means  to  explicitly 
second  kind  expressions  for  the  double  integral  operators  that  arise  when 
analytic  preconditioning  methods  are  applied  to  the  2d  and  3d  scalar  scat¬ 
tering  analogs  of  the  EFIE.  Although  we  ended  up  not  using  these  identities 
in  our  formulation  of  a  well  conditioned  EFIE  for  open  and  closed  surfaces, 
a  brief  discussion  of  the  relevant  identities  is  include  here  for  completeness. 


3.1  2d  scalar 

The  Poincare-Bertrand  identity  is  commonly  stated  as1 


/  dx'  - 1  ^  )  -f  dx"  - 12  —  (; X ")  =  -7T2(/>1  (x)  <j>2  (x)  (x)  (1) 

X  —  XJ_1  x  —  X 

+  -f  dx"(j)2  (x")  l  dx'- - ^  - JJT-Ip  (x") 

J_i  J_i  ( x  —  x')  (ar  —  x")  v  ' 

where  f  means  that  the  integral  is  a  Cauchy  principal  value  integral.  Equa¬ 
tion  (1)  expresses  the  Id  PBI  on  the  open  interval  (—1,1).  A  more  general 

1The  redundant  function  tj)  (x")  is  given  separately  from  <j>2  (x")  for  later  notational 
convenience.  It  will  correspond  to  the  unknown  source  function. 
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form  of  the  Id  PBI  is 


/  **  M  rfO,  (x")  “££>*(**)  _ 


>c 

-  (x)  4*2  (x)  Ip  (x) 


(2) 


+/  dl"4> 2  (x")  -f  dl' (pi  (x') 
Jc  Jc 


r,  dG  (x,  x')  3G  (x',  x") 


dV 


dl" 


ip  (x") 


where  the  integrals  are  now  over  arc  length  on  an  arbitrary  smooth  curve 
C ,  the  derivatives  are  tangential  (to  C)  derivatives,  and  G  is  the  2d  Laplace 
or  Helmholtz  kernel,  i.e., 


G  (x,x) 


-^iogdx-x'l)  or  G  (x, x') 


(3) 


respectively. 

The  2d  analogs  of  the  EFIE  are  a  first  kind  integral  equation 


ipinc  (x)  =  [  dl'G  (x,  x')  ^  for  TM  polarization  (4) 

Jc  “n 

and  a  hypersingular  integral  equation 

=  _  4  r  dl, acix  x-)^  for  TE polarization  (5) 

an  an  J q  on 

where  the  derivatives  are  normal  (to  C)  derivatives. 

If  curve  C  is  closed,  one  obtains  a  second  kind  integral  equation  from  (4) 
or  (5)  by  pre-  or  post-conditioning  it  by  the  integral  operator  of  the  other 
equation  [5]  [6] .  The  resulting  composite  operator  can  be  manipulated  into  a 
sum  of  terms  involving  compact  operators,  plus  a  term  whose  form  matches 
that  of  the  LHS  of  (2)  with  (p±  (x7)  =  (p2  (xw)  =  1.  This  term  has  the  spectral 
characteristics  of  a  second  kind  operator.  Using  (2),  we  can  write  it  in  an 
explicitly  second  kind  form. 

If  C  is  an  open  curve,  one  can  apply  the  same  procedure  to  obtain  second 
kind  integral  equations,  provided  additional  factors  containing  the  appro¬ 
priate  endpoint  singularity  behavior  are  included  in  each  integral  operator 
[7]  [8].  In  the  case  of  the  open  interval  (—1, 1),  this  amounts  to  letting 


(pi  (x')  =  (l  -  (x')2)  /  and  cp2  (x")  =  (l  -  (.'c")2)=F  /  .  (6) 


Since  (pi  (x)  (p2  (x)  =  1  in  both  cases,  the  transformation  of  (2)  again  pro¬ 
duces  an  explicitly  second  kind  integral  operator. 

Notes: 
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•  Having  an  explicitly  second  kind  EFIE  formulation  for  2d  scattering 
is  helpful  because  it  makes  solution  methods  based  on  patch-based 
discretizations  practical,  but  it  is  not  essential  because  in  2d  an  entire- 
domain  discretization  based  on  arc  length  is  often  available. 

•  The  form  of  the  Id  PBI  suggests  a  possible  means  to  obtain  a  second 
kind  integral  equation  for  a  target  with  both  open  and  closed  charac¬ 
teristics.  Consider,  for  example,  a  circle  with  a  thin  fin  attached.  It 
may  be  possible  to  obtain  a  second  kind  formulation  for  this  geometry 
by  choosing  02  (x)  to  be  a  function  that  incorporates  the  appropriate 
endpoint  singularity  behavior  at  the  open  end  of  the  fin  and  smoothly 
approaches  unity  at  the  other  end  of  the  fin,  and  by  choosing  0i  (x) 
to  incorporate  the  reciprocal  behavior. 

•  The  PBI  is  needed  for  computational  purposes  only  when  the  field 
and  source  points  are  in  close  proximity.  Most  interactions  are  far 
interactions  and  they  can  be  computed  using  standard  discretizations 
of  the  integral  operators  or  using  fast  methods  such  as  the  FMM. 


3.2  3d  scalar 


We  can  write  the  2d  extension  to  the  Id  PBI  [2]  [3]  as 

f  ds' 4> l  (x7)  V'G  (x,  x7)  •  f  ds"cj>2  (x,;)  V'G  (x^x77)  0  (x/;) 
Js  Js 


=  (x)  02  (x)0(x)  + 


f  ds"(t> 2  (x7/)  f  ds'(j) i  (x7)  V'G  (x,  x7)  •  V'G  (x7,x77)  0  (x7/)  (7) 

Js  Js 


where  the  integrations  are  over  a  sufficiently  smooth  surface  S  and  X7'G  de¬ 
notes  the  gradient  of  the  3d  Laplace  kernel  with  respect  to  x7.  The  extension 
of  (7)  to  the  case  of  the  3d  Helmholtz  kernel  is  trivial. 

The  equations  for  scalar  scattering  from  a  surface  S  under  Dirichlet  and 
Neumann  boundary  conditions  are  analogous  to  (4)  and  (5),  respectively. 
The  procedures  for  pre-  or  post-conditioning  these  equations  to  obtain  sec¬ 
ond  kind  integral  equations  are  also  analogous  to  those  of  the  2d  scalar 
case. 
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4  Analytically  preconditioning  the  EFIE 

The  electric  field  integral  equation  (EFIE)  for  a  PEC  is 

TJ 

/ - - - * - ; - — n 

— n  (x)xEmc  (x)  =  n  (x)  x  J  ds'  i^ikG  (x,x')  J  (x)  +  j-V  (VG  (x,  x')  •  J  (x')) 

(8) 

where  n  (x)  is  the  unit  normal  at  the  field  point  x,  Emc  (x)  is  the  incident 
electric  field,  k  =  s/JjZuj  is  the  propagation  wavenumber  in  the  interaction 
medium,  G(x,x;)  is  the  Helmholtz  kernel,  and  J  (xr)  is  the  equivalent  sur¬ 
face  source  distribution  (or  current).  The  current  is  related  to  the  surface 
component  of  the  magnetic  field  H  by  J  =  Z n  x  H,  where  Z  =  \fjije  is  the 
wave  impedance  in  the  interaction  medium. 

The  ill  conditioned  nature  of  the  EFIE  is  a  direct  result  of  the  near  field 
behavior  of  the  EFIE  integral  operator,  which  we  can  write  as  the  sum  of 
singular  and  hypersingular  component  operators: 

Ts  =  ikn  (x)  x  f  ds'  G  (x,  x') 

Js 

Th  =  JS(x)xV  /  ds'  VG  (x,x') 

k  Js 

The  singular  component  Ts  is  a  smoothing  operator  -  eigenvalues  of  Ts 
corresponding  to  high  spatial  frequency  eigenmodes  tend  toward  the  origin. 

The  hypersingular  component  TH  is  a  differential  operator  -  eigenvalues  of 
Th  corresponding  to  high  spatial  frequency  eigenmodes  tend  toward  infinity. 

Note  that  Ts  is  the  3d  EM  analog  of  the  TM  operator  on  the  RHS  of  (4) 
and  Th  is  the  3d  EM  analog  of  the  TE  operator  in  (5).  The  EFIE  combines 
both  sources  of  ill  conditioning  into  a  single  integral  equation,  adding  to 
the  challenge  of  finding  a  suitable  analytic  preconditioner  for  open  surface 
targets. 

4.1  Original  analytic  preconditioner 

As  described  in  [1],  the  ideal  preconditioner  for  the  EFIE  for  closed  surface 
targets  is  the  EFIE  operator  itself.  In  other  words,  T 2  =  TT  is  a  second 
kind  operator.  At  first  glance,  this  may  seem  surprising  since  the  expanded 


(9) 

(10) 
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product  operator  contains  a  double  hypersingular  term: 


T2J  =  (Ts  +  Th)  ( Ts 


Th)J 

SrpH 


_  ^  rj-iSrj-iS  _|_  rj-iH  rj-iS  _|_  rj~iS  rpH  _|_  rj~iH  J 


(11) 


However,  integration  by  parts  reveals  that  the  double  hypersingular  term 
actually  vanishes 


ThThJ  =  (n  x  V)  [  ds'  VG  •  (S'  x  V')  [  ds"  V'G  •  J" 

Js  Js 

=  (n  x  V)  [  ds'  GV  ■  (n  x  V')  [  ds"  VG  ■  J"  =  0,  (12) 

Js  ' - . - -  Js 


=  0 


which  leaves 


T2J  =  (ThTs  +  TsTh  +  TSTS)  J  (13) 

The  combination  of  the  first  two  terms  on  the  RHS  of  (13)  is  a  second 
kind  operator.  The  last  term  is  obviously  compact.  With  some  further 
manipulation  (see  [1]),  we  can  re-write  (13)  as 


which  has  the  appealing  property  that  the  kernel  of  each  integral  operator 
involves  no  more  than  one  gradient  on  the  Green  function. 

The  surface  current  J  can  be  decomposed2  into  two  parts,  a  transverse 
part  Jr  that  has  no  surface  divergence  (V  •  J2  =  0)  and  a  longitudinal 
component  JL  that  has  no  surface  curl  ((n  x  V)  •  JL  =  0).  One  can  show 

2 According  to  the  Hclmholtz-Hodge  decomposition  theorem  [9],  a  surface  vector  field 
can  be  decomposed  into  a  curl-free  component,  a  divergence-free  component,  and  a  har¬ 
monic  remainder  (which  is  both  curl-  and  divergence- free).  When  it  exists,  the  harmonic 
component  of  J  needs  to  be  treated  separately. 


analytically  [1]  that  THTS  (or  T°Tt)  is  a  second  kind  operator  for  JT 
and  is  compact  for  JL.  Conversely,  TSTH  (or  T^TL)  is  a  second  kind 
operator  for  JL  and  is  compact  for  Jr.  Therefore  the  sum  THTS  +  TSTH 
(or  TaTT +T@Tl)  is  a  second  kind  operator  for  the  full  current  J  =  + JL. 

4.2  Eigenvalues  and  resonances 

The  magnetic  field  integral  equation  (MFIE), 

n  (x)  x  Hinc  (x)  =  -  J  (x)  -n  (x)  x  V  f  ds'  G  (x,  x')  J  (x')  ,  (15) 

2  Js 

" - V - ' 

I< 

only  applies  to  closed  surface  targets. 

If  S’  is  a  closed  surface,  the  T2  operator  obeys  the  following  identity  from 
Roach  [10]: 

T2  =  -\+K2={-\+K){l+K)-  (16) 

Since  K  is  compact,  the  Roach  identity  confirms  that  T 2  is  a  second  kind 
operator  whose  eigenvalues  accumulate  at  — 

The  factored  representation  of  T 2  in  (16)  shows  that  resonances  of  \  -\~K 
and  —  ^  +  K  are  also  resonances  of  T 2 .  So  even  though  preconditioning  both 
sides  of  (8)  turns  it  into  a  second  kind  integral  equation,  it  does  not,  by 
itself,  always  produce  a  well  conditioned  integral  equation.  There  are  a 
number  of  simple  and  effective  ways  (e.g.,  see  [1]  [1 1] )  to  eliminate  these 
spurious  resonances  without  compromising  the  effectiveness  of  the  analytic 
preconditioner. 

Part  of  the  solution  to  avoiding  spurious  resonances  is  to  use  a  different 
wavenumber  for  the  preconditioning  operator.  For  example,  if  k2  =  k  is 
the  wavenumber  of  the  scattering  problem,  one  can  use  k\  =  ik  for  the 
preconditioning  operator.  As  we  saw  in  our  previous  work  [1],  the  resultant 
integral  operator  T  (k\ )  T  (fo)  is  no  longer  a  single  second  kind  operator, 
but  rather  the  sum  of  two  second  kind  operators  with  different  asymptotic 
collection  point,  one  at  +|  and  another  at  —  In  the  course  of  this  program 
we  realized  that  a  simple  modification  to  the  preconditioning  procedure  puts 
both  asymptotic  collection  points  back  at  —  \ . 

The  sum  of  T  (£2)  preconditioned  by  asTs  (k±)  and  T  (A^)  precondi¬ 
tioned  by  aHTH  (k±)  ,  where  as  and  aH  are  arbitrary  coefficients,  can  be 
rewritten  as 

aH^Ta  (fci)  Tt  (, k2 )  +  ask±TP  (An)  TL  (, k2 )  -  asTs  (An)  Ts  (, k2 ) .  (17) 
k\  k2 
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Setting  as  =  aH  =  1  produces  the  same  result  as  before.  On  the  other 
hand,  setting 

H  k  1  J  S  k 2 

a  =  —  and  a  =  — 
k‘2  k\ 

puts  both  eigenvalue  collections  points  at  —  \ ,  independently  of  the  values  of 
k±  and  .  Having  a  single  asymptotic  eigenvalue  collection  point  improves 
iterative  solver  convergence. 


5  Numerical  validation  of  2d  PBI 


The  2d  Poincare-Bertrand  identity  (PBI)  can  be  written  as 

f  ds'  4>\  (x7)  (n  x  VG  (x,x7))  •  f  ds"  (x77)  (n77  x  V77G  (x7,x77))  ^  (x77) 

Js  Js 

=  |</>i(x)02(x)V’(x)+  (18) 

f  ds"  (j)2  (x77)  if;  (x77)  f  ds1  4>i  (x7)  (n  x  VG  (x, x7))  •  (n77  x  V77G  (x7,x77)) 

Js  Js 


where  S'  is  a  smooth  surface  (open  or  closed),  <f>i ,  <f>2,  and  ^  are  general 
smooth  functions,  and  G  is  the  3d  Laplace  kernel  (i.e. ,  G(x,x7)  =  1/  (47rr), 
where  r  =  |x  —  x7|).  The  identity  expresses  the  fact  that  the  two  integral 
operators 


H\  =  /  ds'  4>\  (x7)  (n  x  VG  (x,x7))  • 

Js 

(19a) 

H2=  [  ds"  <fo  (x77)  (a77  x  VG  (x7,x77)) 

Js 

(19b) 

do  not  commute.  Instead,  the  commutator  is 

[Hi,H2\  =  ^1  (x)  <f>2  (x). 

The  identity  also  holds  for  the  Helmholtz  kernel  (i.e.,  G  (x,  x7)  =  eikr /  (47r r)) 
since  the  Helmholtz  kernel  can  be  decomposed  into  the  sum  of  a  singular 
(Laplace)  part  and  a  regular  part  and  for  the  regular  part  the  operators  do 
commute. 

The  first  objective  of  this  contract  was  to  numerically  demonstrate  the 
validity  of  this  identity  in  the  restricted  case  where  (f>\  =  $2  =  1  and  G  (x,  x7) 
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is  the  3d  Laplace  kernel,  i.e. 


f  ds  (nx  VG  (x,x))'  f  ds"  (n"  x  V"G  (x,x"))  i/i  (x") 

Js  Js 


=  |V’(x)  + 


f  ds"  ^  (x")  f  ds'  (n  x  VG  (x,  x'))  ■  (n"  x  V"G  (x',  x"))  , 

Js  Js 


or 


Idht  (x)  =  (x)  +  Ipbi  V’  (x)  , 

where  Idht  and  Ipsi  represent  the  double  integral  operators,  defined  by 

/dj/t^x)^  [  ds'  (nxVG(x,x))'  f  ds"  (n "  x  V"G  (x  ,x"))  iJj  (x") 

Js  Js 

(20a) 

Ipbi  (x)  =  f  ds"  V’  (xr/)  f  ds1  (n  x  VG  (x,  x’))  ■  (n"  x  V"G  ( x;,x "))  . 
Js  Js 

(20b) 

We  developed  prototype  code  in  C++  that  evaluates  each  of  these  terms 
and  demonstrates  the  validity  of  this  identity  on  several  types  of  surface 
patches. 


5.1  Integral  operator  properties 

We  can  consider  the  validation  procedure  as  involving  the  evaluation  of  four 
different  integrals: 


I\  (x)  =  f  ds1  (nxVG  (x.x'))  1  (x) 

Js 

h  (x)  =  [  ds"  (n"  x  V"G  (x',x"))  /  (x") 

Js 

h  (x)  =  [  ds"  h  (x,  x")  /  (x") 

Js 

I4  ( x,x ')  =  f  ds1  (ii  x  VG  (x,x7))  •  (n”  x  V"G  (x^x")) 

Js 

where  /  (x")  is  an  arbitrary  scalar  function  (which  one  could  later  use  to 
represent  the  unknown  source  distribution)  and  f  (x7)  is  an  arbitrary  surface 
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vector  function.  The  combination  of  I\  and  I2  gives  Idht  and  the  combi¬ 
nation  of  I3  and  1. 4  gives  Ipbi •  The  four  integrals  involve  different  types 
of  kernel  singularities  and  must  be  handled  carefully  to  achieve  accurate, 
numerically  realizable  results. 

Our  prototype  code  evaluates  the  double  Hilbert  transform  integral  Idht  f  (x77) 
by  multiplying  the  separate  matrix  representations  (i.e. ,  discretized  versions) 
of  I\  (x)  and  I2  (x7).  By  contrast,  we  evaluate  the  double  integral  Ipbi  f  (x77) 

“all  at  once”  using  numerical  evaluations  of  the  inner  integral  I4  (x,  x77)  at 
each  quadrature  point  of  the  outer  integral  I3  (x) .  The  “all  at  once”  in¬ 
tegration  procedure  is  more  time  consuming,  but  is  also  more  reliable  for 
computing  accurate  results. 

We  have  observed  the  following  properties: 

•  The  kernels  of  I\  and  I2  have  singularities  that  scale  as  cos  ( 9 )  r-2 
(where  r  is  the  distance  from  the  source  point  to  the  field  point  and 
6  measures  the  angle  around  the  field  point)  so  they  are  integrable, 
but  only  in  the  principal  value  sense.  This  requires  some  care.  For 
example,  one  can  rewrite  each  integral  using  integration  by  parts  as 
the  sum  of  a  edge  integral  of  a  regular  function  and  a  surface  integral 
of  a  function  whose  kernel  singularity  is  no  worse  than  r_1. 

•  ii  (x)  is  a  smooth  function  of  x  for  xsS.  For  x  near  an  open  edge  of 
S,  I\  (x)  varies  as  log  (r)  where  r  is  the  distance  from  x  to  the  open 
edge.  I2  (x7)  behaves  similarly. 

•  The  integrand  of  I4  (x,x77)  is  singular  at  x  and  x77.  Each  singularity 
scales  like  cos  ( 6 )  r-2,  where  r  is  the  distance  from  the  source  point  to 
the  singular  point  and  6  measures  the  angle  around  the  singular  point. 

As  with  I 1  and  I2,  they  are  integrable  singularities  in  the  principal 
value  sense  because  of  the  angular  factor.  When  x  and  x77  are  close 
together  special  care  must  be  taken  to  evaluate  the  integral  accurately 
because  the  integrand  varies  rapidly  between  positive  and  negative 
infinities.  We  have  investigated  several  approaches  to  evaluating  this 
integral,  including  some  that  involve  analytic  singularity  subtraction. 

They  all  work,  but  with  varying  degrees  of  efficiency  and  robustness. 

•  On  a  flat  surface,  I4  (x,  x77)  is  a  smooth  function  of  x77  for  a  given  x, 
even  at  x  =  x".  If  S  is  curved,  then  I4  (x,  x77)  has  a  log  (r)  singularity, 
where  r  =  |x  —  x77|.  For  x77  near  an  open  edge  of  S,  I4  (x,x77)  varies 
as  log  (r)  where  r  is  the  distance  from  x77  to  the  open  edge. 
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•  In  one  sense  /;->  (x)  is  the  easiest  of  the  integrals  to  evaluate  since  the 
singularities  of  its  kernel  I4  (x.  x")  are  no  worse  that  log(r).  How¬ 
ever,  it  can  be  the  most  time  consuming  to  evaluate  because  at  each 
quadrature  point  one  needs  an  accurate  evaluation  of  1 4  (x,  x") . 

5.2  Numerical  results 

Our  prototype  code  can  evaluate  Idht  f  (x)  and  Ipbi  f  (x)  with  adjustable 
accuracy  on  a  variety  of  patch  types  with  a  variety  of  testing  functions.  The 
computation  takes  a  ‘digits’  parameter,  which  specifies  the  accuracy  goal  (in 
number  of  digits)  of  the  adaptive  integration  routines.  For  Idht  f  (x),  there 
is  also  an  ‘order’  parameter.  It  specifies  the  number  of  x.'  points  at  which 
I2  (x')  is  evaluated,  which  is  also  the  number  of  testing  functions  f  (x;)  used 
to  compute  local  corrections  for  I±  (x). 

The  testing  functions  we  chose  are  products  of  polynomials  in  the  surface 
parameterization  u\  and  u2.  The  tables  below  are  organized  such  that  the 
entry  in  the  mth  row  and  nth  column  corresponds  to  the  scalar  function 

/  (x")  =  /  (ui,u2)  =  Pm  ( ui )  Pn  (u2) 

where  Pn  ( u )  is  the  rath  Legendre  polynomial. 

For  easy  comparison,  the  tables  show  computed  values  of  ( Idht  —  \)  f  (x) 
and  Ipbi  /(x). 
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5.2.1  Square  patch 

Tables  1  and  2  give  ( Idht  ~  3)  /  (x)  and  Tables  3  and  4  give  Ipsi  /  (x) 
on  a  square  patch  at  x  =  (0.25,0.35).  Table  5  gives  the  difference  between 
( Idht  ~  \)  f  (x)  (from  Table  2)  and  Ipbi  f  (x)  (from  Table  4). 


Table  1:  ( Idht  — 

I )  /  (*)  for  a 

square  patch  (digits=10  and  order=8) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10514744 

0.02500997 

-0.00878341 

-0.00051174 

-0.01157399 

1 

0.03965702 

-0.01959232 

-0.00205752 

0.00387275 

0.00661861 

2 

-0.02384387 

0.00801955 

0.01640285 

-0.00927404 

0.00050915 

3 

0.00873452 

-0.00266640 

-0.01145919 

0.00844339 

-0.00138054 

4 

-0.01365415 

0.00459554 

0.00706801 

-0.00530567 

0.00024560 

Table  2:  ( Idht  — 

l)  f(x)  for  a 

square  patch  (digits=10  and  order=10) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10517396 

0.02501531 

-0.00882176 

-0.00052061 

-0.01164151 

1 

0.03968532 

-0.01959947 

-0.00199779 

0.00391992 

0.00678708 

2 

-0.02383721 

0.00800491 

0.01638150 

-0.00929691 

0.00045280 

3 

0.00878461 

-0.00267033 

-0.01135429 

0.00839380 

-0.00129506 

4 

-0.01365726 

0.00464442 

0.00699183 

-0.00547012 

0.00012730 

m/n 

Table  3: 
0 

Ipbi  f  (x)  for 
1 

a  square  patch  (digits=6) 

2  3 

4 

0 

-0.10515237 

0.02500882 

-0.00879625 

-0.00052500 

-0.01162961 

1 

0.03967821 

-0.01960292 

-0.00200786 

0.00388514 

0.00675064 

2 

-0.02385828 

0.00801653 

0.01636798 

-0.00930678 

0.00043937 

3 

0.00879925 

-0.00272262 

-0.01136585 

0.00844514 

-0.00131449 

4 

-0.01368469 

0.00462558 

0.00704841 

-0.00545916 

0.00016687 
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Table  4: 

Ipbi  f  (x)  for 

a  square  patch  (digits=8) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10515263 

0.02500873 

-0.00879683 

-0.00052512 

-0.01163036 

1 

0.03967805 

-0.01960298 

-0.00200822 

0.00388503 

0.00675015 

2 

-0.02385886 

0.00801634 

0.01636670 

-0.00930704 

0.00043773 

3 

0.00879902 

-0.00272271 

-0.01136637 

0.00844499 

-0.00131519 

4 

-0.01368545 

0.00462536 

0.00704675 

-0.00545947 

0.00016474 

Table  5:  (. IDHt  )  f  (x)  - 

-  Ipbi  f  ( x ) 

for  a  square  patch 

m/n 

0 

1 

2 

3 

4 

0 

-0.00002133 

0.00000658 

-0.00002493 

0.00000451 

-0.00001115 

1 

0.00000727 

0.00000351 

0.00001043 

0.00003489 

0.00003693 

2 

0.00002165 

-0.00001143 

0.00001480 

0.00001013 

0.00001507 

3 

-0.00001441 

0.00005238 

0.00001208 

-0.00005119 

0.00002013 

4 

0.00002819 

0.00001906 

-0.00005492 

-0.00001065 

-0.00003744 
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5.2.2  Ogive  patch 

Tables  6  and  7  give  ( Idht  ~  3)  /  (x)  and  Tables  8  and  9  give  Ippi  /  (x) 
on  an  ogive  patch  at  x  =  (0.25,  0.35).  Table  10  gives  the  difference  between 
( Idht  ~  f  (x)  (from  Table  7)  and  Ipbi  f  (x)  (from  Table  9). 


Table  6:  ( Idht  ~ 

l)  f(x)  for  an 

ogive  patch  (digits=10  and  order=8) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10497214 

0.01531734 

-0.00434583 

0.00114802 

-0.01055496 

1 

0.02877761 

-0.00922837 

-0.00005219 

0.00023867 

0.00550509 

2 

-0.01299734 

0.00165363 

0.01443126 

-0.00469695 

-0.00057568 

3 

0.00033876 

0.00098081 

-0.00769284 

0.00398050 

-0.00029916 

4 

-0.01035685 

0.00215251 

0.00574442 

-0.00231159 

-0.00099865 

Table  7:  ( Idht  — 

l)  /  0*0  for  an 

ogive  patch  (digits=10  and  order=10) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10496593 

0.01530208 

-0.00434915 

0.00108919 

-0.01056134 

1 

0.02877142 

-0.00924484 

-0.00001177 

0.00045047 

0.00550813 

2 

-0.01300874 

0.00168139 

0.01440599 

-0.00460083 

-0.00046958 

3 

0.00032182 

0.00092830 

-0.00763120 

0.00357058 

-0.00050284 

4 

-0.01027579 

0.00219966 

0.00561910 

-0.00213368 

-0.00077999 

m/n 

Table  8: 
0 

Ipbi  f  (x)  for 
1 

an  ogive  patch  (digits=6) 

2  3 

4 

0 

-0.10469727 

0.01519943 

-0.00457952 

0.00114313 

-0.01047368 

1 

0.02852133 

-0.00909148 

0.00027490 

0.00028032 

0.00536309 

2 

-0.01319548 

0.00178295 

0.01468790 

-0.00468231 

-0.00059081 

3 

0.00047419 

0.00080978 

-0.00791684 

0.00375302 

-0.00029897 

4 

-0.01027494 

0.00217714 

0.00564040 

-0.00215186 

-0.00082333 
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Table  9: 

Ipbi  f  (x)  for 

an  ogive  patch  (digits=8) 

m/n 

0 

1 

2 

3 

4 

0 

-0.10469748 

0.01519944 

-0.00457999 

0.00114316 

-0.01047425 

1 

0.02852126 

-0.00909165 

0.00027474 

0.00028003 

0.00536283 

2 

-0.01319595 

0.00178299 

0.01468690 

-0.00468222 

-0.00059207 

3 

0.00047410 

0.00080952 

-0.00791710 

0.00375258 

-0.00029937 

4 

-0.01027553 

0.00217723 

0.00563913 

-0.00215167 

-0.00082491 

Table  10:  (IDht  ~  |)  f  (x) 

-  IpBI  f  (x) 

for  an  ogive  patch 

m/n 

0 

1 

2 

3 

4 

0 

-0.00026845 

0.00010264 

0.00023084 

-0.00005397 

-0.00008709 

1 

0.00025016 

-0.00015319 

-0.00028651 

0.00017044 

0.00014530 

2 

0.00018721 

-0.00010160 

-0.00028091 

0.00008139 

0.00012249 

3 

-0.00015228 

0.00011878 

0.00028590 

-0.00018200 

-0.00020347 

4 

-0.00000026 

0.00002243 

-0.00002003 

0.00001799 

0.00004492 
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5.3  Discussion 

As  the  ‘digits’  parameter  increases,  the  computed  value  of  Ipbi  /  (x)  con¬ 
verges  for  each  scalar  function.  Likewise,  as  the  ‘order’  parameter  increases, 
the  computed  value  of  Idht  f  (x)  converges  for  each  scalar  function. 

The  computed  values  of  Ipbi  /  (x)  generally  show  faster  convergence 
than  those  of  I dht  f  (x) ,  which  is  a  direct  consequence  of  the  fact  that 
we  are  using  a  two  stage  approach  to  computing  Idht  /(x),  rather  than 
the  all-at-once  approach  used  to  compute  Ipbi  /(x).  If  we  were  to  com¬ 
pute  Idht  /  (x)  using  double  adaptive  integration,  as  we  do  to  compute 
Ipbi  f  (x),  the  convergence  rate  would  be  similar. 

In  any  case,  ( Idht  ~  \)  f  (x)  and  Ipbi  f  (x)  converge  to  approximately 
the  same  value  for  each  patch.  Table  5  shows  that  for  the  square  patch  the 
difference  is  less  than  6  x  10-5  for  all  25  scalar  functions  and  Table  10  shows 
that  for  the  ogive  patch  the  difference  is  less  than  3  x  10~4  for  all  25  scalar 
functions. 

Although  this  does  not  constitute  a  proof  that  the  2d  PBI  identity  is 
correct,  it  does  demonstrate  that  it  is  feasible  to  numerically  evaluate  all  the 
integrals  in  the  identity  and  it  provides  numerical  support  for  the  validity 
of  the  identity. 

6  Analytic  preconditioner  (APEFIE)  for  closed  sur¬ 
face  PEC  targets 

6.1  3d  scalar 

In  the  process  of  numerically  validating  the  2d  PBI  we  came  to  the  realization 
that,  while  the  2d  PBI  affords  the  advantage  of  achieving  an  explicitly  second 
kind  integral  operator,  it  may  still  be  possible  to  construct  a  numerical 
representation  of  the  double  Hilbert  transform  operator  in  such  a  way  that 
its  desirable  spectral  properties  are  not  compromised.  The  advantage  of  the 
latter  approach  is  that  it  is  easier  and  less  time  consuming  to  compute  the 
double  Hilbert  transform  than  the  PBI  integral.  If  we  could  make  do  with 
the  double  Hilbert  transform  alone,  we  could  speed  up  the  calculation  of  the 
near  interaction  matrix  elements  in  the  analytically  preconditioned  EFIE 
considerably. 

Analytically,  the  double  Hilbert  integral  behaves  like  a  second  kind  op¬ 
erator.  Its  spectrum  consists  of  several  distinct  eigenvalues  corresponding 
to  low  frequency  eigenmodes  plus  a  collection  of  eigenvalues  piling  up  at 
The  eigenvalues  at  |  correspond  to  the  higher  frequency  eigenmodes.  The 
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problem  we  previously  observed  was  that  numerical  representations  of  this 
operator  did  not  display  second  kind  behavior.  The  spectrum  of  the  dis¬ 
cretized  double  Hilbert  transform  consisted  of  lower  frequency  eigenvalues 
that  were  well  resolved  and  higher  frequency  ones  that  were  not.  Instead  of 
all  the  eigenvalues  corresponding  to  higher  frequency  eigenmodes  collecting 
at  j,  we  observed  a  trail  of  eigenvalues  extending  from  ^  to  the  origin. 

The  problem  had  to  do  with  the  way  we  were  discretizing  the  double 
Hilbert  transform  term.  The  double  Hilbert  transform  operator  is  the  com¬ 
position  of  two  Hilbert  transforms.  If  we  define  the  individual  Hilbert  trans¬ 
form  operators  as  in  (19),  then  the  double  Hilbert  transform  operator  in 
(20a)  is 

#12  =  H\H2 

In  the  past  we  had  been  computing  the  matrix  representation  of  H\2  (= 
Idht)  as  the  product  of  the  matrix  representations  of  H\  and  H2.  If  S 
represents  a  single  patch  containing  N  sample  points,  then  the  matrix  rep¬ 
resentation  of  H2  was  a  square  matrix  connecting  the  N  source  points  x" 
to  N  intermediate  field  points  x7;  and  H\  was  a  square  matrix  connecting 
the  N  intermediate  source  points  x;  to  the  N  field  points  x.  The  resultant 
product  matrix  H\2  connects  the  N  source  points  x;/  to  the  N  field  points 
x  as  it  should. 

As  part  of  the  process  of  validating  the  2d  PBI,  we  computed  the  matrix 
representation  of  H\2  two  new  ways.  The  first  way  was  to  compute  H\2  ’’all 
at  once”.  In  other  words,  we  computed  the  outer  integral  term  in 

#12  (V’(x))  = 

f  ds'  4>i  (x7)  (n  x  VG  (x, x7))  •  f  ds"  (j)2  (x77)  (n77  x  V77G  (x7,x77))  ip  (x77) 

Js  Js 

(21) 

adaptively  while  also  using  adaptive  integration  to  evaluate  the  inner  integral 
at  each  quadrature  point  of  the  outer  integration.  The  second  way  is  much 
like  the  method  described  in  the  previous  paragraph  except  that  we  used 
M  intermediate  source/field  points  x7  in  the  calculation  of  the  individual 
matrix  representations,  where  M  >  N.  The  matrix  representations  of  H\ 
and  H-2  became  complementary  rectangles  rather  than  squares.  The  final 
result  for  H\2  is,  of  course,  still  square. 

We  computed  the  spectrum  of  H\2  for  a  sphere  using  all  three  tech¬ 
niques.  As  before,  we  found  that  the  original  technique  produced  a  trail  of 
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eigenvalues  from  |  to  the  origin,  which  correspond  to  unresolved  high  spa¬ 
tial  frequency  eigenmodes  on  the  sphere.  When  we  did  the  same  calculation 
using  either  of  the  two  new  ways  this  was  not  the  case.  Instead,  we  saw  the 
correct  spectrum.  For  the  all-at-once  evaluation  method  this  is  as  expected 
-  the  full  operator  is  being  computed  to  whatever  accuracy  is  specified  so 
the  spectrum  should  be  similarly  accurate.  However,  we  found  that  it  also 
worked  for  the  matrix  product  method  when  M  was  only  modestly  larger 
than  N.  In  particular,  we  found  that  doubling  the  sample  point  density  in 
each  surface  dimension  (i.e. ,  M  =  4 N)  was  generally  sufficient  to  recover 
the  correct  spectral  properties.  This  method  of  evaluation  is  generally  much 
faster  than  the  all-at-once  method,  which,  in  turn,  is  much  faster  than  eval¬ 
uating  the  PBI  integral.  When  M  2>  N,  results  obtained  using  the  matrix 
product  method  should  become  indistinguishable  from  those  obtained  using 
the  all-at-once  method  given  that  both  methods  will  be  employing  a  large 
number  of  intermediate  sample  points  distributed  over  the  surface  in  their 
computations. 

The  following  plot  demonstrates  the  difference  between  the  old  and  new 
ways  of  computing  the  double  Hilbert  transform  term.  The  plot  shows 
eigenvalues  for  the  double  Hilbert  transform  operator  on  an  r  =  0.43667A 
sphere.  The  red  points  labeled  “old  way”  were  computed  using  the  matrix 
product  method  with  the  number  of  intermediate  points  equaling  the  num¬ 
ber  of  field/source  points  on  each  patch.  The  green  points  were  computed 
using  the  all-at-once  integration  method.  The  blue  points  are  the  analyti¬ 
cally  computed  eigenvalues3  of  the  double  Hilbert  transform  operator  on  the 
sphere.  The  computation  was  performed  using  a  low  order  discretization  so 
the  red  and  green  eigenvalues  are  seen  to  converge  rather  slowly  to  the  exact 
blue  ones.  The  salient  feature  of  the  plot  is  the  difference  in  the  behavior 
of  the  red  and  green  points  corresponding  to  high  spatial  frequency  modes. 
The  red  eigenvalues  form  a  trail  extending  from  (0.25,0)  to  the  origin.  By 
contrast,  the  green  eigenvalues  tend  to  accumulate  near  (0.25,0),  as  they 
should. 

6.2  3d  vector  EM 

The  3d  vector  EM  analogs  of  the  product  operator  H\H-2  in  (21)  are  the 
product  operators  T°Tt  and  T'T[j  defined  in  (14).  Together  they  constitute 

3  The  Ith  eigenvalue  for  the  double  Hilbert  transform  operator  on  a  sphere  of  radius 
r  is  — l  ( l  +  1)  (Jhfcyfcr))  ,  where  J;  and  H;  are  the  Riccati-Bessel  and  Riccati-Hankel 
functions,  respectively,  of  order  l. 
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Figure  1:  Eigenvalues  spectrum  for  double  Hilbert  transform  operator  on  a 
sphere. 


the  second  kind  part  of  the  T2  operator.  In  our  original  implementation  of 
the  Analytically  Preconditioned  Electric  Field  Integral  Equation  (APEFIE), 
we  obtained  a  matrix  representation  for  T2  by  separately  discretizing  its 
component  operators  Ta,  Tr ,  T'3,  TL,  and  Ts ,  multiplying  them  in  pairs, 
and  adding  them  together.  The  matrix  representation  of  each  operator  was 
square.  The  numerical  spectrum  exhibited  eigenvalues  tending  toward  the 
origin  just  as  we  saw  in  the  scalar  case  because  the  discretized  operator  was 
under-resolving  high  spatial  frequency  eigenmodes  that  the  discretization 
was  capable  of  representing. 

In  our  original  implementation  we  overcame  this  problem  by  adding  the 
magnetic  field  equation  (MFIE)  to  the  TTIE4  to  create  a  combined  field 
equation  (CFIE).  This  not  only  moves  the  low-pass  filtered  eigenvalues  of 
T 2  away  from  the  origin  by  virtue  of  the  constant  term  (one  half)  contributed 

4TTIE  is  an  alternative  abbreviation  for  APEFIE.  The  name  evokes  the  fact  that  it  is 
based  on  the  TT  operator. 
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by  the  MFIE,  it  also  eliminates  spurious  resonances.  In  the  present  case,  we 
were  able  to  solve  the  problem  by  slightly  over-resolving  the  discretization  at 
intermediate  points  of  the  product  operator  like  we  did  for  the  double  Hilbert 
transform  operator  in  the  scalar  case.  This  approach  has  the  advantage  that 
it  can  also  be  extended  to  the  open  surface  case  where  adding  the  MFIE  is 
not  an  option. 

Figures  2,  3,  and  4  compare  the  spectrum  of  the  TTIE  computed  the  old 
way  (with  the  number  of  intermediate  sample  points  equal  to  the  number 
of  field/source  points)  and  the  new  way  (with  the  four  times  a  many  inter¬ 
mediate  sample  points  as  field/source  points).  In  these  plots,  the  red  points 
are  for  the  old  discretization  method  (a  standard  iterative  EFIE  solver), 
the  green  points  are  for  the  new  discretization  method  (more  intermediate 
sample  points),  and  the  blue  points  are  the  analytical  result.  Of  the  three 
targets  one  expects  the  spectral  behavior  of  the  APEFIE  to  be  ideal  only 
for  the  sphere  because  its  surface  is  smooth.  By  contrast,  the  cube  has 
sharp  edges  and  corners  and  the  ogive  has  pointy  tips,  all  of  which  are  ge¬ 
ometric  singularities.  The  ideal  APEFIE  for  targets  with  such  geometric 
singularities  probably  incorporates  complementary  singularity  functions  in 
the  APEFIE  operator  (see  discussion  in  the  next  section).  Although  we  did 
not  do  that,  the  computed  spectra  of  the  APEFIE  operators  for  the  cube 
and  ogive  have  the  characteristics  one  desires  in  the  ideal  operator,  namely 
a  bounded  spectrum  and  an  eigenvalue  collection  point  well  separated  from 
the  origin. 

7  Analytic  preconditioner  for  open  surface  PEC 
targets 

The  APEFIE  overcomes  the  local  causes  of  EFIE  ill  conditioning  on  smooth 
targets  for  which  the  source  distribution  is  also  expected  to  be  smooth  and 
continuous.  Open  surfaces  are  not  smooth  targets  and  one  expects  the  source 
distribution  to  exhibit  singular  behavior  near  an  open  edge.  Therefore,  it 
should  not  be  surprising  if  the  ideal  analytic  preconditioner  for  open  surface 
PEC  targets  turns  out  to  be  different  from  the  APEFIE  for  closed  surface 
PEC  targets.  In  fact,  the  quest  for  the  ideal  APEFIE  for  open  surfaces  may 
provide  some  insight  into  what  the  ideal  APEFIE  for  closed  surface  targets 
with  geometric  singularities  (edges,  corners,  sharp  tips,  etc.)  should  look 
like. 

To  get  some  insight  into  what  the  ideal  APEFIE  for  open  surface  targets 
should  look  like,  it  is  helpful  to  first  analyze  the  problem  of  scalar  scattering 
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PEC  sphere  «i^«nv^|jee 
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Figure  2:  Eigenvalue  spectrum  for  the  APEFIE  on  a  PEC  sphere  with  radius 
r  =  j A,  comparing  old  (red)  and  new  (green)  discretization  methods. 


from  open  surface  targets  in  2d  and  3d.  In  both  cases,  analytical  solutions 
are  available  for  scattering  from  simple  canonical  targets.  The  solutions  are 
particularly  simple  in  the  Laplace  limit,  which  is  sufficient  for  analyzing  the 
local  behavior  of  the  operators.  Unfortunately,  there  is,  as  yet,  no  simple 
analytical  solution  for  the  3d  EM  vector  case5. 

7.1  Analytical  solutions 

7.1.1  2d  scalar  case 

The  2d  scalar  analogs  of  the  EFIE  are  (4)  and  (5)  and  the  canonical  open 
surface  target  for  2d  scalar  scattering  is  a  flat  strip. 

’There  are  series  solutions  (e.g.,  [12])  for  scattering  from  a  flat  PEC  disc  but  they  are 
quite  complicated  and  do  not  offer  immediate  insight  into  the  modes  of  the  EFIE  operator 
on  the  disc. 
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Figure  3:  Eigenvalue  spectrum  for  the  APEFIE  on  a  cube  with  1A  sides, 
comparing  old  (red)  and  new  (green)  discretization  methods. 


Dirichlet  boundary  conditions 

First  consider  (4),  which  is  the  appropriate  integral  equation  for  2 d  scalar 
scattering  from  the  curve  C  under  Dirichlet  boundary  conditions  (ip  =  0  on 
C )  or,  equivalently,  TM-polarized  3d  vector  EM  scattering  from  the  infinite 
cylinder  with  cross  section  C.  Preconditioning  (4)  by  the  hypersingular 
kernel  integral  operator  from  (5)  followed  by  an  arbitrary  scalar  function 
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PEC  ogive  eigenvalues 
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Figure  4:  Eigenvalue  spectrum  for  the  APEFIE  on  a  PEC  ogive  with  tip-to- 
tip  length  l  =  10A  and  center  radius  r  =  1A,  comparing  old  (red)  and  new 
(green)  discretization  methods. 
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-  dl  jcdl  —ai^h  (x)  Jcdl  G(x-x )  -^r 


where  we  have  also  made  use  of  the  2d  identity  [13] 

dd G(x,x7)  2  ,  n  9  3G(x,x') 

- - =  (a-a)k°(x’x)-gi—^r- 


dn  dn' 


(22) 


(23) 
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Using  integration  by  parts  we  can  replace  the  last  term  on  the  far  RHS  of 
(22)  by 


_d_  [  u,dG  (x,x;) 

dl  Jc  dl 


=  /  dl 

Jc 


+  /  dl 

Jc 


91  jc 

,d G  (x,x7)  d0 i  (x') 

dl  Jl 

d G  (x,  x') 

~ai 


7— -01  (x')  j  dl"G  (x,  x") 

-01  (x)  J  dl 


dip  (x") 

,dG(x,x'),  ,  A  f  „<9G(x',x")  d0(x") 


c 


dl' 

dl"G  (x',x") 


dn" 

dip  (x") 
dn" 


(pi  (x;)  J  dl"G  (x',x") 


dip  (x") 
dn" 


x'=C2 


(24) 


J  x'=Ci 


Let  us  specify  C  to  be  the  interval  from  —1  to  1  and  let  the  Green  func¬ 
tion  be  the  2d  Laplace  kernel  given  in  (3).  For  later  notational  convenience, 
we  will  also  pull  an  arbitrary  scalar  function  02  ( x ")  out  of  he.,  re¬ 

place  d^n  ^  by  02  (x")  £  ( x ").  Then  the  first  term  on  the  far  RHS  of  (24) 
becomes 


1 

W) 


-  x'J- 1 


dx" 


02 {x") 

x'  —  x" 


t  A") 


(25) 


which,  apart  from  the  constant  prefactor,  is  identical  to  the  LHS  of  (1). 

At  this  point  it  is  worthwhile  to  make  note  of  the  following  standard 
identities  ((22.13.4)  and  (22.13.3)  in  [14]): 


dx  V7!  x'2Un- 1  (x  )  -  7rTn  (x) 

x  —  x'  v  ' 

(26) 

dx  - (x  )  =  TlUn  1  (x) 

X  -  x'  VI  -  X/2 

(27) 

where  Tn  and  Un  represent  the  nth  Chebyshev  polynomials  of  the  first  and 
second  kind,  respectively,  and  n  =  1,2,.... 

If  we  define 


01  (x) 


yj  1  —  x2  and  02  (x)  =  1/ \/l  —  x2 


(28) 


it  follows  that 


1 

(2^? 


01  (aQ  f 
~  x'J_  1 


dxw 


02  (•?") 

x'  —  x" 


Tn  (x") 


(29) 
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for  n  =  1,2,3,....  In  other  words,  with  this  choice  for  <f>i  and  fa,  the 
eigenfunctions  of  the  operator  in  (25)  are  the  Chebyshev  polynomials  of  the 
first  kind  and  in  all  cases  the  eigenvalue  is  —1/4. 

Notes: 


•  Any  smooth  function  defined  on  the  interval  [—1  :  1]  can  be  expressed 
as  a  linear  combination  of  first  kind  Chebyshev  polynomials; 


•  The  product  of  4> 2  (%)  with  any  linear  combination  of  first  kind  Cheby¬ 
shev  polynomials  will  have  l/\/l  =p  x  singularities  at  the  endpoints  ±1 
and  be  regular  in  the  interior  of  the  interval,  which  matches  the  end¬ 
point  behavior  one  expects  for  solutions  to  (5); 

•  4>\  (x)  vanishes  at  the  endpoints. 


The  identity  for  XJn  ( x )  corresponding  to  (29), 


(2tt)2/-i  x-x'J_  i 


dx>' ~7~~JJUn  (x")  =-\Un  (*)  >  (30) 


applies  for  n  =  0,1,2 
singularity  functions. 


Note  the  reversed  order  of  the  endpoint 


•  Equations  (29)  and  (30)  are  special  cases  of  (1)  for  which  the  last  term 
on  the  RHS  of  (1)  evaluates  to  zero. 


Since  there  is  a  smooth  mapping  from  the  interval  [—1  :  1]  to  any  smooth 
curve  C,  we  can  deduce  that  the  double  integral  operator 


IC 


jj/  dG  (x,  x')  A 
di - ^ - (t> i  (x  ) 


/  c 


j,n  d°  (x',  x")  (  , 

dl  - Qfl - ^2  (*  ) 


(31) 


is  a  second  kind  operator.  All  the  other  terms  in  our  operator  decomposi¬ 
tion  are  compact.  Since  the  difference  between  the  Laplace  kernel  and  the 
Helmholtz  kernel  is  a  regular  function,  this  result  also  holds  for  the  integral 
operators  associated  with  2d  scattering.  Therefore,  in  the  general  case,  one 
could  write  a  second  kind  integration  equation  for  open  surface  Dirichlet 
scattering  as 


d_ 

dn 


dl' 


ic 


dn 


ic 


dl"G  (x',x")</> 2  (x")£(x"), 


(32) 
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where  (pi  and  <p2  are  complementary  endpoint  singularity  functions  like  those 
in  (28)  and  the  solution  is  given  by  ^  =  4>2  (x")  £  (x,r)-  The  eigenvalues 
for  high  spatial  frequency  eigenmodes  of  the  double  integral  operator  on  the 
RHS  of  (32)  will  accumulate  at  —1/4. 

Note: 

•  If  the  curve  C  was  the  unit  circle  (a  closed  surface)  instead  of  the 
unit  interval  (an  open  surface),  a  similar  analysis  results  in  the  same 
conclusion  for  Eq.  (32).  The  only  difference  is  that  the  functions  (p\ 
and  cp2  both  equal  one.  The  same  analysis  is  extensible  to  the  case 
where  C  is  any  smooth  closed  curve. 


Neumann  boundary  conditions 

A  similar  result  obtains  for  2d  scalar  scattering  from  the  curve  C  under  Neu¬ 
mann  boundary  conditions  ( dnip  =  0  on  C)  or,  equivalently,  TE-polarized 
3d  vector  EM  scattering  from  the  infinite  cylinder  with  cross  section  C.  In 
this  case,  however,  it  is  more  convenient  to  make  use  of  a  postcondition¬ 
ing  operator  that  contains  the  complementary  kernel  singularity  and  the 
complementary  singularity  factor6. 

It  works  like  this.  Assume  that  we  can  write  ^  as 

ip  (x')  =  <p\  (x')  /  dl"G  (x',  x")  (p2  (x")  £  (x")  (33) 

Jc 

for  some  functions  (pi  (x')  and  £  (x/7),  in  which  case  we  can  rewrite  (5)  as 


dipmc  (x) 
dn 


d_  f  ,dG  (x,xQ 
dn  Jc  dn' 


<Pi  (x)  ^  dl"G  (x,  x")  (p2  (x")  £  (x")  . 


(34) 


After  applying  (23)  and  using  integration  by  parts  as  before,  we  can 

6  Preconditioning  by  the  same  operator  would  also  result  in  a  second  kind  operator  for 
ip.  However,  one  cannot  manipulate  it  into  a  sum  of  compact  operators  and  a  double 
Hilbert  transform  operator  because  the  integration  by  parts  procedure  produces  a  term 

whose  integrand  is  proportional  to  d</)(/  ^  ~  (l  —  x ,2)  3,/2,  which  is  not  an  integrable 
endpoint  singularity. 
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rewrite  the  RHS  of  (34)  as 


/c  M  Jc  <r  M  { (x») 


+ /c  j-8G(87x')a^r')  x  d,"c  <x'  x">  *  <x"> « <x") 

—  k2  f  dl '  (n  •  n7)  G  (x,  x7)  0i  (x7)  f  dl" G  (x7,  x77)  02  (x77)  £  (x77) 

7C  Jc 

'dG  (x,  x7 


-01  (x7)  /  dZ77G  (x,x")  02  (x")  £  (x") 

Jc  J  x'=Ci 


x'=C2 


(35) 


Following  the  same  arguments  as  in  the  Dirichlet  case  we  find  that  the 
first  term  in  (35)  is  second  kind  with  respect  to  ^  whereas  the  other  terms 
are  compact  operators  or  zero.  Therefore  (34)  is  a  second  kind  integral 
equation  for  open  surface  Neumann  scattering.  Having  obtained  a  solution 
for  £,  the  solution  for  0  follows  immediately  from  (33). 


7.1.2  3d  scalar  case 

The  equations  describing  scalar  wave  scattering  under  Dirichlet  and  Neu¬ 
mann  boundary  conditions  in  3d  are  the  same  as  (4)  and  (5)  except  that  the 
integrals  over  curve  C  must  be  replaced  by  integrals  over  surface  S.  The 
canonical  open  surface  target  for  3d  scattering  is  a  flat  circular  disc. 

The  analytical  preconditioning  procedure  is  similar  to  what  we  did  in 
the  2 d  scalar  case,  i.e.,  pre-  or  postcondition  each  integral  equation  by  an 
integral  operator  that  includes  the  complementary  kernel  singularity  and  the 
edge  singularity  factor  that  is  complementary  to  the  expected  source  term 
edge  singularity.  The  rest  of  the  procedure  follows  much  the  same  path  as 
before.  The  important  differences  are  in  the  details  of  the  second  kind  term, 
which  we  can  analyze  in  the  Laplace  limit  for  the  specific  case  where  S  is 
the  unit  disc. 

Oblate  spheroidal  coordinates  are  ideal  for  analyzing  the  disc  problem 
because  the  3d  scalar  Helmholtz  and  Laplace  equations  are  separable  in  this 
coordinate  system  {£,  77,  ip }  and  the  degenerate  (£  — ►  0)  limit  of  an  oblate 
spheroid  is  a  flat  disc.  Two  views  of  the  oblate  coordinate  system  on  a  disc 
are  shown  in  Fig.  5.  The  two  oblate  spheroidal  coordinates  that  describe  a 
disc  are  cp  and  r).  The  p  coordinate  is  the  usual  azimuthal  coordinate.  On  a 
disc  of  radius  a,  the  7]  coordinate  is  related  to  the  radius  by  r  =  a  (l  —  r/2) . 

Note: 
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Side  view 


Top  view 


Figure  5:  Oblate  spheroidal  coordinate  system 


(  2\1/2 

•  Since  77  =  ( 1  —  (r/a)  )  ,  one  expects  solutions  to  the  Dirichlet  prob¬ 

lem  to  scale  like  1  / 77  near  the  edge  of  the  disc  and  solutions  to  the  Neu¬ 
mann  problem  to  scale  like  77  near  the  edge  of  the  disc.  In  other  words, 
77  and  1  / 77  can  serve  as  complementary  edge  singularity  functions. 

•  As  77  runs  from  —1  to  0  to  +1,  the  location  on  the  disc  runs  from  the 
underside  center,  to  the  edge,  to  the  topside  center. 

Solutions  to  the  Laplace  equation  on  the  disc  are  linear  combinations  of 
modes  of  the  form 

(36) 

where  P™  (77)  is  the  associated  Legendre  function  of  order  ( l,m ).  The  func¬ 
tion  P,m  (77)  is  an  even  (odd)  polynomial  in  77  times  rm  for  l  +  m  even  (odd). 

Solutions  to  the  Helmholtz  equation  with  wavenumber  k  on  the  disc  are 
linear  combinations  of  modes  of  the  form 

'film  (77,  <P]  7)  =  S™  (17, 7)  eirmf>  (37) 

where  S'”1  (77, 7)  is  the  oblate  spheroidal  wavefunctions  of  order  ( l,m )  and 
dimensionless  scale  size  7  =  ka.  Unlike  the  associated  Legendre  functions, 
the  oblate  spheroidal  wavefunctions  do  not  have  a  simple  closed  form.  Nor 
do  they  obey  most  of  the  simple  identities  that  the  associated  Legendre 
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polynomials  obey  (see  8.73-8.74  of  [15]).  In  the  low  frequency  limit,  they 
reduce  to  associated  Legendre  functions,  i.e. , 

as  7-0. 

The  key  ingredients  for  the  analysis  are: 

•  Angular  orthogonality: 


2t r 


j dip  =  2tt 5mim, 


Radial  orthogonality  ((7.112.1)  of  [15]): 

>M  /  \  DH  r*A  _  2  +  lmD' 


(??)  pr  w  = 


-1 


2/  +  1  (l  —  |m|)! 


Su 


Differential  area  element: 


dA  =  a2  j??]  dr]  dip 

2? r  1 


(Note  that  a2  J  dip  J  \rj\  dr]  covers  the  disc  twice  -  once  over  the  top 


0  -1 

side  and  once  over  the  bottom  side.) 


Laplace  Green  function  expansion  (derived  from  (10.3.63)  of  [16])  in 
terms  of  Laplace  modes  ^p|m  (77)  eimv^j : 


a  (*p = t  <■«>  pr  w 


47raf-^  ^  cim  (Z+M)! 

1=0  m=— l 
(by  2) 


(38) 


where  a  is  the  radius  of  the  disc  and 


C-lm  — 


r  (7+^+2)  r  (i^^n+2) 


(39) 
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With  these  relations  we  obtained  the  following  analytical  results7  when 
S  is  disc  of  radius  a  and  G  is  the  Laplace  kernel: 


•  First  kind  and  hypersingular  integral  operators  including  edge  singu¬ 
larity  factors: 


a 


V’Zm  (V,  <P) 


for  l+m  even  (40) 


and 


d  f  A  /  dG  (x,  x^)  , 


V  rf  ) 


Clm  fpl  m 
2  a  rj 


for  l  +  m  odd 


(41) 


where  ipim{7h(f)  is  the  scalar  Laplace  mode  as  defined  in  (36).  Note 
that  i/jim  (rj,  ip)  is  an  even  (odd)  function  of  rj  if  l  +  m  is  even  (odd),  so 
tplm  (rj,  ip)  in  (40)  with  l  +  m  is  even,  and  ^’lrndh(p)  jn  ^4^  with  l  +  m 
is  odd,  are  both  regular  functions  on  the  whole  of  the  disc  (i.e.,  for 
—  1  <  77  <  1  and  0  <  <p  <  27 r)  and  even  functions  of  rj. 

•  Combinations  of  first  kind  and  hypersingular  integral  operators  includ¬ 
ing  edge  singularity  factors: 

d  f  .  fdG  (x,x  )//',//„//  //\  1  ,  /  n  rr\ 

l 

=  X  alm^l'm  Ob  +)  for  1  +  m  even  (42) 

l'=m 
(by  2) 


and 


-  /  ds'  G  (x  x')  f  ^ lm  y//) 

Js  ’  rj'  dm'  Js  dn"  \  rj" 

=  X  /?L#m(7?,y)  for  Z  +  to  odd 

l'=m+ 1  ^ 

(by  2) 

where  afm  and  are  (unspecified)  constant  coefficients. 


(43) 


'  As  far  as  we  know,  these  identities  are  novel  results. 
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The  integral  operator  in  (42)  does  not  mix  m' s.  For  a  given  (/,  m)  with 
l  +  m  even,  the  eigenfunctions  of  this  operator  are  linear  combinations 
of  ipi'm.  {Vi  v)  with  l'  =  l,  l  —  2, , . . .  m  .  The  eigenvalue  corresponding 
to  the  ( l,m)th  eigenfunction  is 

1(1  +  l)2  -  m2 


n2 

“7m 


(44) 


Likewise,  the  integral  operator  in  (43)  does  not  mix  m’s.  For  a  given 
( l,m )  with  l  +  m  odd,  the  eigenfunctions  of  this  operator  are  linear 
combinations  of  - with  l ’  =  l,  l  —  2, , . . .  m  +  1.  The  eigenvalue 
corresponding  to  the  ( l,m)th  eigenfunction  is 


1  Clm 

4  l2  —  m2 


(45) 


Combinations  of  Hilbert-transform-like  integral  operators  including 
edge  singularity  factors: 

f  ds  n  x  VG  (x,  xr)  •  r\  f  ds"  n;  x  V'G  (x7,  x")  —i/iim  {jf ,  +") 

Js  Js  V 

1 1  (l  + 1)  -  m2  ,  .  .  . 

iplm  (■ rj ,  99)  tor  l  +  m  even  (4b) 

i’lm  W,  <P" 


r2 

~'lm 


and 


f  ds'  SxVG  (x,  x')  •  *  f  ds "  S'  x  V'G  (x',  x")  r,"  f 
J s  V  Js  V 

1  {l(l  +  !)  -  ™?)  CL  V’ lm  (V,  <p) 


V 


^  ((l  +  l)2  —  ra2^  ( l 2  —  m 2)  d 


for  l  +  m  odd, 


(47) 


Dirichlet  boundary  conditions 


The  integral  equation  appropriate  to  Dirichlet  boundary  conditions  is  (4) 
with  curve  C  replaced  by  surface  S.  Preconditioning  both  sides  of  this 
equation  by  the  surface  form  of  the  hypersingular  integral  operator  from  (5) 
followed  by  an  arbitrary  surface  function  <f>i  (xr)  gives  the  surface  analog  to 
(22),  namely 


+  Ss  ds"G  *x',x”^2  (X"K(X"). 


(48) 
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where  we  have  substituted  <f>2  fa")  £  (x")  for  d ^  .  Letting  fa  (xr)  =  rf 
and  fa  fa")  =  1  /rf  puts  the  RHS  of  (48)  into  the  same  form  as  the  integral 
operator  in  (42).  The  eigenfunctions  of  this  operator  are  regular  functions 
of  position  on  the  disc.  The  added  factor  of  fa  fa")  =  1  /rf  gives  them 
the  correct  edge  singularity  behavior  for  scalar  solutions  to  the  Laplace  or 
Helmholtz  problem  on  the  disc  under  Dirichlet  boundary  conditions.  The 
eigenvalues  of  the  integral  operator  in  (48)  in  the  Laplace  limit  can  be  de¬ 
duced  from  (44). 

Alternatively,  one  can  take  (48)  a  step  further,  using  the  identity  [13] 


d  dGfa,x') 
dn  dn' 


(n  •  n')  k2G  (x,  x')  —  (ii  x  V)  •  (n;  x  V7)  G  (x,  x')  (49) 


and  performing  an  integration  by  parts  as  was  done  in  the  2d  scalar  case. 
The  resulting  RHS  is  the  sum  of  a  compact  operator  and  the  term8 


(n  x  VG  (x.x7))  -fa  (x7)  j  ds"  (n;  x  V'G  (x^x"))  fa  (x7/)  £  (x7/) 


where  fa  fa'),  fa  fa"),  and  £  (x/;)  are  defined  as  before.  The  eigenfunctions 
and  corresponding  eigenvalues  of  this  operator  in  the  Laplace  limit  are  given 
by  (46). 


Neumann  boundary  conditions 

A  similar  result  obtains  for  Neumann  boundary  conditions.  There  are  two 
ways  to  go  about  it.  One  can  either  precondition  the  integral  equation  ap¬ 
propriate  to  Neumann  boundary  conditions  by  the  integral  operator  from 
the  Dirichlet  problem  or  use  the  post-conditioner  approach  described  previ¬ 
ously  for  2d  scalar  scattering  under  Neumann  boundary  conditions. 

In  the  first  case,  the  integral  operator  consists  of  a  compact  term  and  a 
term  of  the  form 

-  js  ds'  G  (x’  x)  ^2  (x)  ^  Js  ds"dG(^n;rX  Vi  (x")  £  (x")  (51) 

where  i\)  (x")  has  been  replaced  by  fa  fa")  £  (x").  Letting  fa  fa")  =  if  and 
fa  fa')  =  1/?/  puts  (51)  into  the  same  form  as  the  operator  in  (42),  whose 
eigenvalues  (in  the  Laplace  limit)  are  given  by  (45).  In  the  second  case,  the 

8This  operator  in  differs  from  the  double  Hilbert  transform  operator  —H12  defined  in 
(21)  by  at  most  a  compact  operator. 
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integral  operator  consists  of  a  compact  term  and  a  term  identical  to  that 
of  (50),  where  (x7)  and  <f>2  (x")  are  defined  as  above.  In  both  cases,  the 
added  factor  of  (j)\  (x")  =  rf  in  the  expression  for  the  solution  guarantees 
that  it  exhibits  the  correct  edge  singularity  behavior  for  the  solutions  under 
Neumann  boundary  conditions. 

Second  kind  character  of  the  eigenvalue  spectrum 


l 


Figure  6:  Spectrum  of  the  3d  scalar  double  integral  operator  in  (42)  in  the 
Laplace  limit;  i.e.,  a  plot  of  the  eigenvalues  given  in  (44). 


The  eigenvalue  spectrum  in  the  3d  scalar  case  is  not  as  simple  as  it  was  in 
the  2d  scalar  case  wherein  the  eigenvalue  associated  with  every  eigenfunction 
was  —  In  the  3d  scalar  case,  it  depends  on  l  and  m.  For  example,  Figure 
6  is  a  plot  of  the  Dirichlet  eigenvalues9  in  the  Laplace  limit  for  0  <  m  <  l 

9  By  virtue  of  the  identity  the  eigenvalue  spectrum  for  the 

Neumann  case  as  a  function  of  l  +  1  is  the  same  as  the  Dirichlet  spectrum  as  a  function  of 
L,  i.e.  the  Neumann  counterpart  to  Figure  6  would  look  the  same  except  that  all  points 
would  be  shifted  rightward  by  one  unit.  The  spectrum  of  (50)  is  the  same  as  the  Dirchlet 
spectrum  shown  in  Figure  6. 
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up  to  l  =  50.  In  the  high  spatial  frequency  limit  (l  — »  oo),  the  eigenvalues 
asymptotically  approach  collection  points  whose  values  depend  on  l  —  m. 
The  number  of  collection  points  is  infinite,  but  they  all  fall  into  the  range 
from  —  f  —  —0.393  to  —  Obviously,  this  does  not  conform  to  the  classic 
definition  of  a  second  kind  spectrum  because  the  eigenvalues  do  not  tend  to 
accumulate  at  a  single  point.  It  is  more  like  a  sum  of  second  kind  operators 
all  of  whose  eigenvalue  collection  points  are  between  —  and  The 

important  point  for  our  purposes  is  that  the  spectrum  is  bounded  and  all 
the  eigenvalues  stay  well  away  from  the  origin  in  the  high  spatial  frequency 
limit. 

The  corresponding  spectrum  for  the  Helmholtz  case  will  be  different  from 
that  shown  above,  especially  for  l  <  ka.  However,  the  high  spatial  frequency 
( l  — >  oo)  behavior  of  the  modes  in  the  Helmholtz  case  is  the  same  as  that 
in  the  Laplace  case  so  the  above  conclusions  regarding  the  limiting  behavior 
of  the  eigenvalues  and  the  second  kind  nature  of  the  operator  apply  equally 
well  to  3d  scalar  scattering. 

7.1.3  3d  vector  EM  case 

Preconditioning  the  3d  EM  vector  case  is  more  complicated  than  either  of 
the  previous  scalar  cases  because  the  EFIE  operator  has  both  Dirichlet  and 
Neumann  characteristics.  The  Ts  component  of  T  behaves  like  the  Dirichlet 
operator  in  (4)  and  the  TH  component  of  T  behaves  like  the  Neumann 
operator  in  (5).  The  vector  source  distribution  also  contains  both  Dirichlet 
and  Neumann  characteristics.  Near  an  open  edge,  the  component  of  the 
current  parallel  to  the  edge  obeys  Jj  oc  d-1/2  like  the  Neumann  source 
distribution  and  the  component  of  the  current  perpendicular  to  the  edge 
obeys  Jj_  oc  5 1/2  like  the  Dirichlet  source  distribution,  where  5  is  the  distance 
to  the  edge. 

We  know  how  to  analytically  precondition  the  integral  equations  for  3d 
scalar  scattering  under  Neumann  or  Dirichlet  boundary  conditions,  so  we 
might  expect  some  combination  of  these  approaches  should  work  in  the  3d 
EM  vector  case.  However,  it  is  not  obvious  how  to  do  this  given  the  way  the 
EFIE  ties  them  together.  Even  so,  we  can  make  some  educated  guesses. 

Informed  by  the  Id  and  2d  PBI  and  the  ideal  analytic  preconditioners 
for  2d  and  3d  scalar  open  surface  scattering,  we  expect  that  the  ideal  an¬ 
alytic  preconditioner  for  3d  vector  EM  open  surface  scattering  should  look 
something  like  this 

“T2” J  =  (Ts  +  Th )  fa  (Ts  +  Th )  02J  (52) 
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where  (f>i  and  <f> 2  are  complementary  edge  singularity  functions  [and  is 
the  source  distribution  including  edge  singularity  behavior]. 

We  evaluated  two  versions  of  “T2”  J  for  open  surfaces  -  one  in  which  <f>  1  = 
1  and  another  in  which  cf>  1  contained  the  complementary  edge  singularity  to 
(j> 2.  In  both  cases,  we  enforced  the  correct  edge  singularity  behavior  on 
the  currents  by  incorporating  it  into  the  discretization.  For  a  Nystrom 
discretization  such  as  ours,  this  is  done  by  locating  sample  points  according 
to  a  quadrature  rule  that  is  high  order  for  integrating  the  appropriate  set 
of  singular  functions.  For  example,  on  quadrilateral  patches  with  one  open 
edge,  we  can  create  a  suitable  2d  quadrature  rule  from  the  tensor  product 
of  Id  quadrature  rules,  one  of  which  is  a  regular  rule  (e.g.,  Gauss-Legendre) 
and  the  other  is  a  rule  for  integrating  polynomials  times  a  weighting  function 
with  a  e>-1/2  singularity  at  one  end  of  the  interval  (e.g.,  Gauss- Jacobi).  One 
can  generate  other  tensor  product  rules  for  cases  where  more  than  one  edge 
touches  an  open  surface10.  We  have  used  these  rules  for  the  past  several 
years  in  our  general  purpose  RCS  code  (FastScat™)  to  achieve  high  order 
convergence  on  open  surface  scattering  problems  with  the  EFIE. 

We  evaluated  the  properties  of  the  APEFIE  operator  in  (52)  on  several 
simple  open  surfaces,  two  of  which  are  discussed  below.  In  our  numerical 
experiments  the  two  versions  of  “T2”  J  behaved  similarly. 

Triangle-circle  target 

The  triangle-circle  (or  EMCC  wedge-cylinder-plate)  target  is  composed  of 
an  equilateral  triangle  attached  to  a  semicircle.  Our  seven  quadrilateral, 
handmade  mesh11  is  shown  in  Fig.  7.  The  discretization  of  each  outer 
quadrilateral  patches  is  based  on  a  product  rule  with  one  singular  edge;  the 
discretization  for  the  center  quadrilateral  is  based  on  a  tensor  product  of 
regular  quadratures.  The  total  (one-sided)  surface  area  is  3.31A2. 

Figs.  8  and  9  show  the  eigenvalue  spectra  for  the  open  surface  specific 
APEFIE  given  in  (52)  and  the  standard  EFIE  in  (8),  for  independent  dis¬ 
cretization  densities  ranging  from  20.6  unknowns/A  to  26.7  unknowns/A.  A 

1()When  two  adjacent  edges  of  a  quadrilateral  touch  open  edges  the  actual  geometry 
includes  a  corner.  We  do  not  know  of  a  high  order  expansion  for  the  currents  at  a 
corner  and  our  product  rule  quadratures  are  unlikely  to  be  ideal.  However,  our  experience 
shows  that  a  tensor  product  of  5±x'2  quadratures  performs  much  better  than  low  order 
quadrature  rules. 

11  This  mesh  exemplifies  one  of  the  differences  between  our  Nystrom  discretization  and 
a  typical  Galerkin  discretization.  Being  patch  based  rather  than  edge  based,  a  Nystrom 
discretization  allows  for  a  vertex  in  the  middle  of  an  edge. 


37 


Figure  7:  Triangle-circle  mesh. 


blue  dot  on  top  of  a  green  dot  on  top  of  a  red  dot  indicates  an  eigenvalue  of 
a  fully  resolved  eigenmode. 

Several  features  are  noteworthy.  The  eigenvalues  of  the  APEFIE  are 
confined  to  a  relatively  small  region  on  the  negative  real  half  of  the  com¬ 
plex  plane  and  tend  to  cluster  about  a  few  points.  With  few  exceptions12, 
the  eigenvalues  have  converged.  The  smallest  eigenvalues  are  still  well  sep¬ 
arated  from  the  origin.  By  contrast,  the  EFIE  spectrum  is  spread  over  a 
wider  region  centered  about  the  origin  on  both  halves  of  the  complex  plane 
and  is  expanding  with  discretization  refinement.  The  APEFIE  spectrum  is 
consistent  with  a  well  conditioned  operator,  the  EFIE  spectrum  is  not. 

This  interpretation  is  borne  out  by  the  iterative  solver  convergence  re¬ 
sults  shown  in  Figs.  10  and  11  and  summarized  in  Table  11.  The  first 
plot  shows  residual  error  vs.  iteration  count  for  the  transpose- free  quasi¬ 
minimum  residual  (TFQMR)  iterative  solver.  The  second  plot  shows  the 
same  information  for  the  multiple- RHS  generalized  minimum  residual  (MGCR) 
solver.  Note  that  axes  scaling  is  log  -  log.  In  all  cases  we  attempted  to  solve 
for  the  RCS  in  both  polarizations  over  the  azimuthal  range  from  0°  to  180° 
in  1°  steps  at  10°  elevation.  A  block  diagonal  preconditioner  based  on  a 
block  size  of  was  used  with  the  EFIE  to  accelerate  iterative  solver 
convergence. 

1JThe  curved  tail  of  unresolved  eigenvalues  may  be  connected  with  the  fact  that,  al¬ 
though  our  discretization  is  high  order  for  integrating  currents  in  the  interior  of  the  patch 
and  near  smooth  portions  of  the  edge,  it  is  not  high  order  near  corners. 
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Figure  8: 


Open  surface  APEFIE  spectrum  for  triangle-circle  target. 


The  differences  in  iterative  solver  convergence  for  the  analytically  precon¬ 
ditioned  EFIE  and  the  block  diagonally  preconditioned  EFIE  are  dramatic. 
The  APEFIE  consistently  converged  to  a  residual  error  of  10~5  in  about  30 
iterations  with  the  TFQMR  and  CGS  solvers.  The  EFIE  obtained  a  solu¬ 
tion  for  the  first  set  of  excitations  in  about  50K  iterations  and  was  unable 
to  reach  the  termination  condition  on  the  second  pass. 

The  EFIE  fared  somewhat  better  with  the  MGCR  solver  and  the  APE¬ 
FIE  fared  somewhat  worse,  but  the  APEFIE  was  still  almost  a  factor  of  10 
faster  to  converge. 


TTIE 

EFIE 

TFQMR 

~30 

~60000 

CGS 

~30 

~50000 

MGCR 

~200 

~1500 

Table  11:  Number  of  iterations  to  convergence. 
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Figure  9:  Standard  EFIE  spectrum  for  triangle  circle  target. 
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RCS  results  from  the  APEFIE  and  EFIE  for  an  H-polarized  incident  field 
are  compared  to  the  EMCC  measurement  results  in  Figs.  12  and  13.  Both 
methods  are  converging  to  the  same  results,  although  the  APEFIE  results 
are  slightly  more  accurate  for  any  given  discretization  density.  All  RCS 
results  should  be  symmetric  about  0°  since  the  target  is  symmetric  about 
0°.  The  EMCC  measurement  data  is  not  symmetric,  however,  indicating 
some  form  of  measurement  error.  We  have  plotted  the  two  halves  of  the 
measurement  data  on  the  same  0°  —  180°  scale  to  provide  a  crude  indication 
of  measurement  accuracy. 

Our  implementation  of  the  Nystrom  method  does  not  enforce  current 
continuity  between  patches.  Instead,  current  continuity  is  a  natural  con¬ 
sequence  of  achieving  high  order  convergence  to  a  exact  solution  that  is 
continuous  on  a  smooth  surface.  This  is  generally  easier  to  achieve  with  a 
second  kind  integral  operator  such  as  the  MFIE  and  than  with  an  integro- 
differential  operator  such  as  the  EFIE. 

Current  amplitudes  on  the  circle-triangle  target  computed  using  the 
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Figure  10:  Iterative  solver  convergence  on  triangle-circle  target  using 
TFQMR  solver  and  EFIE  (with  a  block  diagonal  preconditioner)  and  the 
open  surface  APEFIE. 


APEFIE  and  EFIE  are  shown  in  Figs.  14  and  15.  The  discretization  is 
the  same  in  both  cases:  1694  unknowns  (basis  order  11).  Both  current  dis¬ 
plays  exhibit  the  expected  wave  pattern  and  current  singularities  near  edges. 
The  main  difference  is  that  the  boundaries  between  patches  that  are  easily 
visible  in  the  EFIE  case  are  virtually  invisible  when  the  APEFIE  is  used. 

A  Nystrom  discretization  computes  the  current  only  at  selected  sample 
points  and  we  determine  the  currents  elsewhere  by  fitting  the  sample  point 
data  to  a  high  order  interpolation  function  based  on  the  underlying  quadra¬ 
ture  rule  for  each  patch.  Currents  at  a  boundary  between  patches  may 
be  extrapolated  from  either  side.  The  superior  continuity  of  the  APEFIE- 
computed  currents  at  artificial  boundaries  is  a  consequence  their  Nystrom- 
sampled  currents  being  more  accurate,  which,  in  turn,  is  a  consequence  of 
the  APEFIE  behaving  like  a  second  kind  operator. 
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Figure  11:  Iterative  solver  convergence  on  triangle-circle  target  using  MGCR 
solver  and  EFIE  (with  a  block  diagonal  preconditioner)  and  the  open  surface 


APEFIE. 


Circular  flat  PEC  disc  target 

By  virtue  of  its  symmetry  and  simplicity,  the  circular  flat  PEC  disc  is  the 
canonical  open  surface  PEC  target.  We  investigated  the  properties  of  the 
APEFIE  for  discs  of  different  sizes.  In  general,  the  spectral  characteristics 
of  the  APEFIE  on  a  disc  are  similar  to  those  of  the  APEFIE  on  the  triangle- 
circle  target,  in  that  the  eigenvalues  are  confined  to  a  relatively  small  region 
on  the  negative  real  half  of  the  complex  plane  and  they  tend  to  cluster  about 
a  few  points. 

We  were  surprised,  however,  to  find  that  a  few  of  the  eigenvalues  converge 
to  the  origin  when  the  radius  r  obeys  the  condition  that  kr  is  a  root  of  Jn , 
the  Bessel  function  of  integer  order  n.  The  multiplicity  of  the  eigenvalues 
at  the  origin  is  one  for  n  =  0  and  two  otherwise.  For  example,  the  second 
zero  of  J\  is  7.0156.  The  spectrum  of  the  APEFIE  for  a  disc  of  radius 
R  =  7.0156/27T,  is  shown  in  Fig.  16.  The  dot  at  the  origin  represents  two 
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Figure  12:  Triangle-circle  target  RCS:  EMCC  measurement  vs.  APEFIE 
for  three  different  discretizations. 


converged  eigenvalues. 

How  do  we  interpret  these  zero  eigenvalues?  Do  they  correspond  to 
physical  resonances  on  the  disc?  That  possibility  is  contradicted  by  the  ob¬ 
servation  that  the  resonance  eigenvalues  are  not  apparent  when  we  use  the 
standard  EFIE.  Are  they  spurious  resonances?  It  is  well  known  that  solu¬ 
tions  to  the  EFIE  for  closed  targets  are  susceptible  to  contamination  from 
spurious  resonances,  which  correspond  to  solutions  to  the  interior  scattering 
problem.  Such  solutions  are  also  eigenmodes  of  the  exterior  EFIE  operator 
with  zero  eigenvalues.  They  do  not  contribute  to  the  exterior  field.  However, 
an  open  surface  such  as  the  disc  does  not  have  an  interior.  The  observed 
resonances  seem  likely  to  be  spurious  resonances  of  the  two  open  surface 
APEFIEs  we  have  investigated,  although  we  do  not  have  a  simple  explana¬ 
tion  for  why  they  exist  or  why  they  appear  under  the  particularly  simple 
conditions  described  above. 

In  the  hope  of  identifying  the  underlying  cause  of  the  resonance  eigen¬ 
values  and  perhaps  modifying  our  APEFIE  formulation  to  avoid  them,  we 
attempted  to  find  an  analytical  series  solution  to  the  PEC  disc  scattering 
problem.  We  sought  a  series  solution  (analogous  to  the  Mie  series  solution 
for  scattering  from  a  sphere)  in  terms  of  eigenmodes  of  the  APEFIE  opera- 
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Figure  13:  Triangle-circle  target  RCS:  EMCC  measurement  vs.  EFIE  for 
three  different  discretizations. 


tor.  We  generated  a  basis  set  of  vector  modes  on  the  disc  in  terms  of  oblate 
spheroidal  wavefunctions,  namely 

Xim  (' rj ,  p  ])  =  nx  V  ipim  (v,  P  l) 

Vim  Ob  P  7)  =  V  (V,  P  l)) 

where  ^/m  (rj,  ip\  7)  is  the  3d  scalar  Helmholtz  mode  defined  in  (37)  and  l  +  m 
is  odd.  The  xim  vector  basis  functions  can  represent  surface  vector  fields 
with  no  surface  divergence  whose  azimuthal  (radial)  component  scales  like 
r]-1  ( rj )  near  the  edge.  The  vim  vector  basis  functions  can  represent  surface 
vector  fields  with  no  surface  curl  whose  radial  (azimuthal)  component  scales 
like  r]  (rj3)  near  the  edge.  The  3d  vector  Helmholtz  equation  is  not  separable 
in  spheroidal  coordinates,  so  it  is  not  obvious  what  the  best  choice  of  vector 
basis  set  is. 

We  derived  a  set  of  conditions  for  the  series  solution  coefficients  in  terms 
of  some  generic  Id  integrals  of  products  of  oblate  spheroidal  wavefunctions 
that  one  can  evaluate  numerically  for  any  given  value  of  the  dimensionless 
scale  size  7. 

Late  in  the  program,  we  realized  that  it  would  be  possible  to  analyti- 
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Figure  14:  Current  distribution  on  triangle-circle  target  computed  using 
APEFIE. 


cally  determine  whether  or  not  our  APEFIE  was  a  second  kind  operator  by 
analyzing  its  behavior  in  the  k  — ►  0  limit  where  the  oblate  spheroidal  wave- 
functions  reduce  to  analytically-expressible  associated  Legendre  functions. 
The  keys  to  achieving  second  kind  behavior  are  the  terms 

-  f  ds'  nx  VG  (x,  x;)  rj  f  ds"  n'x V'G  (x7,  x") -n"x  V"ipim  {rf ,  7)  (53) 

V  Js  Js 

and 

V  nx  f  ds'  n'x  VG  (x,x7)  f  ds”  VG  (x',  x")  •  V"  (77",  ip"\ 7))  . 

Js  V  Js 

(54) 

Using  Mathematica,  we  were  able  to  analytically  evaluate  expressions 
(53)  and  (54)  for 

lMr/V';7  =  0  )=Pr(r,)eim* 

with  l  +  m  odd.  We  found  that  for  m  =  0  ,  the  spectra  of  both  operators 
are  second  kind  with  an  asymptotic  eigenvalue  collection  point  at  —  For 
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Figure  15:  Current  distribution  on  triangle-circle  target  computed  using 
EFIE. 


all  other  values  of  m,  however,  the  range  of  each  operator  was  bigger  than 
the  domain,  so  it  was  not  a  closed  system  and  we  were  not  able  to  construct 
eigenfunctions  with  our  basis  set. 

Nonetheless,  we  did  investigate  a  modification  to  the  APEFIE  that 
showed  some  promise  for  eliminating  resonance  eigenvalues.  Not  only  must 
the  tangential  component  of  the  electric  field  vanish  on  the  surface  of  a 
PEC,  so  also  must  the  normal  component  of  the  magnetic  field.  Therefore, 
the  source  distribution  must  obey  the  following  constraint  on  the  normal 
component  of  the  magnetic  field: 

a .  Hinc  =  -a  •  nscat  =  -a  •  [  ds'  vg  (x,  x')  x  j  (x') .  (55) 

Js 

At  any  given  source  point,  the  APEFIE  applies  two  constraints  on  the  two 
components  (one  for  each  surface  direction)  of  the  source  distribution.  Eq. 
(55)  places  an  additional  constraint  on  the  two  source  components.  One 
could  solve  the  over-constrained  system  of  equations  described  by  the  APE¬ 
FIE  and  (55)  (i.e. ,  three  constraints,  two  unknowns)  or  solve  the  evenly 
determined  system  of  equations  formed  by  adding  (55)  to  one  or  both  of  the 
APEFIE  component  equations  (i.e.,  two  constraints,  two  unknowns).  We 
investigated  the  latter  approach. 
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Figure  16:  Eigenvalues  spectrum  for  APEFIE  on  a  PEC  disc  with  a  “reso¬ 
nant”  radius. 


The  spectrum  of  the  resultant  equation  is  shown  in  Fig.  17.  Adding 
the  normal  magnetic  field  constraint  has  the  salutary  effect  of  shifting  the 
resonance  eigenvalues  away  from  the  origin.  However,  it  also  results  in  a 
more  widely  dispersed  spectrum,  which  is  likely  to  have  a  deleterious  effect 
on  the  iterative  converge  rate.  Nonetheless,  this  version  of  the  APEFIE 
represents  a  major  improvement  over  the  conventional  EFIE  in  terms  of 
condition  number  and  iterative  solver  reliability. 

8  Discussion 

Although  we  seem  not  to  have  found  the  ideal  analytic  preconditioner  for  the 
EFIE  on  open  surface  PEC  targets,  we  have  found  one  that  is  practical  to 
implement  numerically  and  that  drastically  improves  the  condition  number 
of  the  linear  system.  A  few  of  the  practical  implications  and  applications  of 
this  achievement  are  discussed  below. 
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Figure  17:  Eigenvalue  spectrum  of  APEFIE  +  constraint  on  the  normal 
component  of  the  magnetic  field. 


8.1  Implications 

•  Discretization  independence.  Although  our  implementation  is  based  on 
a  Nystrom  discretization,  it  could  equally  well  have  been  done  using  a 
different  discretization,  such  as  Galerkin. 

•  Time  savings.  The  purpose  of  this  program  is  to  investigate  ana¬ 
lytic  preconditioners  for  open  surface  PEC  targets  and  demonstrate 
on  small  test  targets  how  an  APEFIE  improves  iterative  solver  per¬ 
formance  as  a  consequence  of  regulating  the  operator  (i.e. ,  improving 
its  condition  number).  The  real  domain  of  practical  applicability  of 
an  APEFIE  is  on  large  targets  for  which  iterative  solvers  in  conjunc¬ 
tion  with  fast  methods  are  a  practical  necessity  and  the  time  spent  in 
the  iterative  solver  phase  dominates  the  overall  solution  time.  Many 
targets  of  practical  interest  fall  into  this  category. 

•  Near  field  accuracy  improvements.  Currents  computed  using  the  EFIE 
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usually  show  the  most  visible  discontinuities  at  patch  boundaries  even 
under  conditions  where  the  far  field  has  converged  to  high  accuracy. 
The  underlying  reason  is  that  the  EFIE  operator  is  hypersingular. 
The  APEFIE,  being  based  on  a  second  kind  operator,  generally  pro¬ 
duces  noticeably  smoother  surface  currents.  Smoother  currents  should 
improve  the  accuracy  of  derived  near  field  quantities  such  as  port 
impedances. 

•  Dual  surface  method  applied  to  the  EFIE.  The  dual  surface  integral 
equation  method  [11]  is  a  straightforward,  general  technique  for  elim¬ 
inating  spurious  resonances  in  closed  surface  scattering  problems.  It 
does  so  by  adding  a  term  derived  from  a  particular  compact  operator 
to  the  original  equation.  There  was  no  particular  benefit  in  using  this 
method  with  the  EFIE  because  the  resulting  equation  would  still  be  ill 
conditioned.  However,  it  does  make  sense  with  the  APEFIE  because 
the  APEFIE  is  well  conditioned  and  it  stays  that  way  even  after  the 
dual  surface  component  is  added  on.  Thus  a  dual  surface  APEFIE 
is  worth  considering  as  an  alternative  to  the  CFIE  on  closed  surface 
PEC  targets. 

8.2  Applications 

•  EM  scattering  calculations.  The  APEFIE  is  useful  for  solving  scat¬ 
tering  problems  involving  both  open  and  closed  PEC  targets,  but  is 
especially  valuable  on  large  open  surface  targets  due  to  the  lack  of 
practical  alternatives. 

•  Quasi-static  EM  analysis.  Since  the  APEFIE  is  well  conditioned  from 
high  frequencies  down  to  DC,  it  could  be  used  to  analyze  quasi-static 
EM  problems  (such  as  motors)  for  which  the  EFIE  is  not  a  good  anal¬ 
ysis  tool  because  of  ill  conditioning. 

•  Scalar/ acoustic  scattering.  It  is  obvious  from  the  analysis  of  the  flat 
circular  disc  in  Section  7.1.2  how  to  construct  an  ideal  analytic  pre¬ 
conditioner  for  scalar  scattering  from  an  open  surface  target  under 
Dirichlet  or  Neumann  boundary  conditions.  One  could  solve  large, 
open  surface  acoustic  scattering  problems  involving  acoustically  hard 
or  soft  surfaces  this  way. 
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8.3  Conclusions 


A  long-standing  problem  with  the  integral  equation  approach  to  modeling 
very  thin  metal  structures  in  the  frequency  domain  is  that  the  choice  of 
integral  equation  is  limited  to  the  inherently  ill  conditioned  EFIE.  For  large 
scattering  problems,  one  would  like  to  be  able  to  speed  up  the  calculation 
by  using  iterative  solvers  in  conjunction  with  fast  solver  methods  (such  as 
the  FMM),  but  the  requirement  to  use  the  EFIE  makes  this  impractical  in 
many  cases.  The  general  purpose  matrix  preconditioners  that  have  been 
applied  to  address  this  problem  operate  with  varying  degrees  of  reliability 
and  computational  efficiency. 

The  ideal  solution  would  be  an  inherently  well  conditioned,  second  kind 
integral  equation  suitable  for  use  on  open  surface  PEC  targets.  The  analytic 
preconditioner  developed  and  investigated  under  this  program  comes  very 
close  to  attaining  this  goal.  Depending  on  your  point  of  view,  we  have 
either  found  a  well  tuned  preconditioning  operator  for  the  EFIE  or  a  new, 
well  conditioned  integral  equation  to  augment  the  EFIE. 

We  implemented  an  analytically  preconditioned  EFIE  (APEFIE)  that 
stabilized  and  dramatically  improved  iterative  solver  performance  as  com¬ 
pared  to  conventional  numerical  preconditioning  methods  when  applied  to 
simple  test  targets.  For  example,  the  number  of  iterations  required  to 
achieve  solution  convergence  on  the  EMCC  triangle-circle  target  fell  from 
~50000  for  the  standard  EFIE  with  a  conventional  block  diagonal  precon¬ 
ditioner  to  ~30  with  the  APEFIE.  Similar  results  were  obtained  on  other 
small  open  surface  test  targets.  The  most  significant  benefits  in  terms  of 
computational  savings  will  be  observed  when  an  fast  method  version  of  the 
APEFIE  is  applied  to  large  open  surface  targets. 

The  underlying  reason  for  the  iterative  solver  performance  improvement 
can  be  found  in  the  eigenvalue  spectrum  of  the  integral  operator.  The  APE¬ 
FIE  has  a  second  kind  spectrum,  with  eigenvalues  corresponding  to  high 
spatial  frequency  eigenmodes  tending  to  cluster  about  a  few  collection  points 
well  separated  from  the  origin.  At  the  origination  of  this  program  we  sus¬ 
pected  that  it  would  be  a  nontrivial  task  to  actually  achieve  such  a  spectrum 
in  a  numerical  implementation  and  we  proposed  to  use  a  2d  version  of  the 
Poincare-Bertrand  identity  (PBI)  to  make  this  task  easier.  Early  on  we  nu¬ 
merically  validated  the  2d  PBI,  but  as  the  program  progressed,  we  found 
an  alternative,  less  cumbersome  means  to  achieve  the  same  objective  so  we 
abandoned  the  2d  PBI  approach. 

Our  investigation  showed  that  under  some  circumstances  the  APEFIE 
spectrum  can  include  a  small  number  of  zero  eigenvalues,  which  are  reminis- 
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cent  of  the  spurious  resonances  known  to  plague  the  EFIE  for  closed  surface 
PEC  targets.  To  understand  this  problem  we  attempted  to  obtain  an  an¬ 
alytical  solution  for  the  spectrum  of  the  APEFIE  on  the  canonical  open 
surface  target,  a  circular  disc.  We  made  some  progress  in  this  direction  (as 
reported  here),  but  in  the  end,  did  not  fully  realize  this  goal  or  arrive  at  a 
satisfactory  understanding  of  the  reason  for  the  zeros  in  the  spectrum.  One 
hopes  that  a  more  complete  analysis  will  bring  forth  an  even  better  APEFIE 
for  which  these  zero  eigenvalues  do  not  exist. 
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List  of  Acronyms 


EFIE 

MFIE 

APEFIE 

TTIE 

FMM 

EMCC 

PBI 

RCS 

PEC 


Electric  Field  Integral  Equation 
Magnetic  Field  Integral  Equation 

Analytically  Preconditioned  Electric  Field  Integral  Equation 
T-squared  Integral  Equation  (synonymous  to  APEFIE) 

Fast  Multipole  Method 
Electromagnetic  Code  Consortium 
Pointcare-Bertrand  Identity 
Radar  Cross-Section 
Perfect  Electrical  Conductor 
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